October 8, 2015

Learned today a bit about the way Astrobear generates a forcing function to perturb the grid in turbulence simulations. The grid is transformed into Fourier space so that each cell in Fourier space is given a wave vector. The wave vectors (with corresponding amplitude and phase information) will be summed to give a forcing function. The amplitudes are semi random in that they are weighted against the spectral index of the power spectrum one wishes to achieve (the amplitude is always orthogonal to the k vector, but how it is orthogonal is random). The phases are random. Once this spectrum is generated in Fourier space it is transformed back into physical space and saved into a variable for the forcing function. This forcing function (really an acceleration function) is then multiplied to every zone of the box for a specified fluid variable (by default the velocity). The forcing function is normalized so that it it produces purely solenoidal perturbations (i.e. the amplitude vector is perpendicular to the k vector). At each time step the velocity is shifted so that the average velocity is 0.

The forcing function goes something like:

\vec{a}(x,y,z) = \sum_{k}{A_{\vec{k}}e^{i(\vec{k} \cdot \vec{x})}}

where

\vec{A}_{\vec{k}} is a random complex vector subject to

|\vec{A}_{\vec{k}}| \propto k^{\beta} for k < k_{\max}

and

\vec{A}_{\vec{k}} \cdot \vec{k} = 0

We also want \vec{a} to be real, so

\vec{A}_{\vec{-k}} = \vec{A}_{\vec{k}}

However, in the code - it looks like both are chosen randomly, and then the real part of the transform is taken. I have to think about the implications of this…

\vec{a}(x,y,z) = \mbox{REAL} \left ( \sum_{k}{\vec{A}_{\vec{k}}e^{i(\vec{k} \cdot \vec{x})}} \right )

\vec{a}(x,y,z) = \sum_{k}{ \mbox{REAL} \left ( \vec{A}_{\vec{k}}e^{i(\vec{k} \cdot \vec{x})} \right ) }

a^j(x,y,z) = \sum_{k}{ |A^j{\vec{k}}|\cos{(\vec{k} \cdot \vec{x}+\theta^j_{\vec{k}})} }

where A^j_{\vec{k}} = |A^j_{\vec{k}}| e^{-i \theta^j_{\vec{k}}}

Now when we combine the \pm \vec{k} terms, we get

|A^j_{\vec{k}}| e^{-i \theta^j_{\vec{k}}} + |A^j_{\vec{-k}}| e^{-i \theta^j_{\vec{-k}}}