While reviewing the source code recently to trace the historical context of the algorithm of Turbsim, I came across a specific scaling factor change in the Coh2H subroutine made about 12 years ago. I am hoping to get some clarification on whether this was an intentional physical modification, or if there might be another reason for this change.
Theoretical Context
In the Veers wind field simulation method, the Fourier coefficients at a specific frequency are constructed as:
$$
V = H \Phi
$$
Where:
-
$H$ is the lower triangular matrix derived from the Cholesky factorization.
-
$\Phi$ is the diagonal matrix of random phases ($e^{i\phi_j}$).
Because the Power Spectral Density (PSD) matrix represents the expected energy (statistical expectation), the cross-spectral matrix $\Sigma$ is theoretically defined by the expected value of the Fourier coefficients:
$$
E[V V^H] = E[H \Phi \Phi^H H^T]
$$
Since the random phases $\phi_j$ are independent and uniformly distributed, $E[\Phi \Phi^H] = I$ (the identity matrix). This is mathematically equivalent to evaluating the matrix structure assuming the phase expectation is unity (i.e., normalizing $\phi = 0$ for the energy verification). Thus, the core requirement is:
$$
H H^T = \Sigma
$$
Code Implementation
The target cross-spectral matrix $\Sigma$ is defined element-wise as $\Sigma_{ij} = \sqrt{s_i s_j} \cdot C_{ij}$, where $s_i$ is the auto-spectrum and $C$ is the coherence matrix. The code first performs Cholesky factorization on $C$ to get the lower triangular matrix $L$ (stored in TRH).
Then, the following code block scales the elements of $L$ to construct the final matrix $H$:
DO J = max(1,p%usr%NPoints),p%grid%NPoints ! Column
TRH(Indx) = TRH(Indx) * SQRT( ABS( S(IFreq,J,IVec) ) ) ! (A) Diagonal elements (I == J)
Indx = Indx + 1
DO I = J+1,p%grid%NPoints ! Row
! S(IFreq,I,IVec) should never be less than zero, but the ABS makes sure...
!Indx = NPoints*(J-1) - J*(J-1)/2 + I !Index of H(I,J)
! TRH(Indx) = TRH(Indx) * SQRT( ABS( S(IFreq,I,IVec) ) ) ! (B) old code 12 years ago
TRH(Indx) = TRH(Indx) * SQRT( SQRT( ABS( S(IFreq,I,IVec) * S(IFreq,J,IVec) ) ) ) ! (C) Off-diagonal elements (I > J)
Indx = Indx + 1
ENDDO !I
ENDDO !JLooking at the scaling factors applied here:
- (A) For diagonal elements ($i = j$): The scaling factor is $\sqrt{s_j}$. This perfectly matches the theoretical requirement $H_{jj} = \sqrt{s_j} \cdot l_{jj}$.
- (B) Old code for off-diagonal elements ($i > j$): The scaling factor was $\sqrt{s_i}$.
- (C) Current code for off-diagonal elements ($i > j$): The scaling factor is $(s_i s_j)^{1/4}$.
Mathematical Derivation ($2 \times 2$ example)
To verify if the current off-diagonal scaling factor $(s_i s_j)^{1/4}$ satisfies $H H^T = \Sigma$, let's look at a $2 \times 2$ matrix.
Let the coherence matrix be $C$ and its Cholesky factorization be $L$:
$$
C = L L^T =
\begin{bmatrix} l_{11} & 0 \\
l_{21} & l_{22}
\end{bmatrix}
\begin{bmatrix} l_{11} & l_{21} \\
0 & l_{22}
\end{bmatrix}
$$
Let $D^{1/2}$ be the diagonal matrix of auto-spectra:
$$
D^{1/2} =
\begin{bmatrix} \sqrt{s_1} & 0 \\
0 & \sqrt{s_2}
\end{bmatrix}
$$
The expected cross-spectrum is $\Sigma = D^{1/2} C D^{1/2} = (D^{1/2} L) (D^{1/2} L)^T$.
Therefore, the theoretical target matrix $H_{theory}$ is:
$$
H_{theory} = D^{1/2} L =
\begin{bmatrix}
\sqrt{s_1} l_{11} & 0 \\
\sqrt{s_2} l_{21} & \sqrt{s_2} l_{22}
\end{bmatrix}
$$
However, using the current code's scaling factor for the off-diagonal element, the constructed matrix $H_{current}$ becomes:
$$
H_{current} =
\begin{bmatrix} \sqrt{s_1} l_{11} & 0 \\
(s_1 s_2)^{1/4} l_{21} & \sqrt{s_2} l_{22}
\end{bmatrix}
$$
Let's check the off-diagonal element (Row 2, Col 1) of the resulting matrix $H_{current} H_{current}^T$:
$$
(H_{current} H_{current}^T)_{21} = H_{current, 21} \cdot H_{current, 11}
$$
$$
= \left( (s_1 s_2)^{1/4} l_{21} \right) \cdot \left( \sqrt{s_1} l_{11} \right)
$$
Since $l_{11} = 1$ and $l_{21} = c_{21}$:
$$
= c_{21} \cdot s_1^{3/4} s_2^{1/4}
$$
But the physically expected off-diagonal element $\Sigma_{21}$ should be:
$$
\Sigma_{21} = \sqrt{s_2 s_1} \cdot c_{21} = c_{21} \cdot s_1^{1/2} s_2^{1/2}
$$
It appears that $s_1^{3/4} s_2^{1/4} = s_1^{1/2} s_2^{1/2}$ only holds true when $s_1 = s_2$. For spectra with $s_1 \neq s_2$, the current off-diagonal scaling factor seems to introduce a slight deviation in the cross-spectrum.
I am curious if this change 12 years ago was made to address a specific numerical stability issue or a different interpretation of coherence back then.
While reviewing the source code recently to trace the historical context of the algorithm of Turbsim, I came across a specific scaling factor change in the
Coh2Hsubroutine made about 12 years ago. I am hoping to get some clarification on whether this was an intentional physical modification, or if there might be another reason for this change.Theoretical Context
In the Veers wind field simulation method, the Fourier coefficients at a specific frequency are constructed as:
Where:
Because the Power Spectral Density (PSD) matrix represents the expected energy (statistical expectation), the cross-spectral matrix$\Sigma$ is theoretically defined by the expected value of the Fourier coefficients:
Since the random phases$\phi_j$ are independent and uniformly distributed, $E[\Phi \Phi^H] = I$ (the identity matrix). This is mathematically equivalent to evaluating the matrix structure assuming the phase expectation is unity (i.e., normalizing $\phi = 0$ for the energy verification). Thus, the core requirement is:
Code Implementation
The target cross-spectral matrix$\Sigma$ is defined element-wise as $\Sigma_{ij} = \sqrt{s_i s_j} \cdot C_{ij}$ , where $s_i$ is the auto-spectrum and $C$ is the coherence matrix. The code first performs Cholesky factorization on $C$ to get the lower triangular matrix $L$ (stored in
TRH).Then, the following code block scales the elements of$L$ to construct the final matrix $H$ :
Looking at the scaling factors applied here:
Mathematical Derivation ($2 \times 2$ example)
To verify if the current off-diagonal scaling factor$(s_i s_j)^{1/4}$ satisfies $H H^T = \Sigma$ , let's look at a $2 \times 2$ matrix.
Let the coherence matrix be$C$ and its Cholesky factorization be $L$ :
Let$D^{1/2}$ be the diagonal matrix of auto-spectra:
The expected cross-spectrum is$\Sigma = D^{1/2} C D^{1/2} = (D^{1/2} L) (D^{1/2} L)^T$ .$H_{theory}$ is:
Therefore, the theoretical target matrix
However, using the current code's scaling factor for the off-diagonal element, the constructed matrix$H_{current}$ becomes:
Let's check the off-diagonal element (Row 2, Col 1) of the resulting matrix$H_{current} H_{current}^T$ :
Since$l_{11} = 1$ and $l_{21} = c_{21}$ :
But the physically expected off-diagonal element$\Sigma_{21}$ should be:
It appears that$s_1^{3/4} s_2^{1/4} = s_1^{1/2} s_2^{1/2}$ only holds true when $s_1 = s_2$ . For spectra with $s_1 \neq s_2$ , the current off-diagonal scaling factor seems to introduce a slight deviation in the cross-spectrum.
I am curious if this change 12 years ago was made to address a specific numerical stability issue or a different interpretation of coherence back then.