Commit e4e2cd61 authored by Pierre Augier's avatar Pierre Augier
Browse files

Better test spatiotemp-spectra-2d

parent 9441d558460d
......@@ -688,24 +688,32 @@ class SpatioTemporalSpectraNS:
# get kz, kh
params_oper = self.sim.params.oper
deltakz = 2 * pi / params_oper.Lz
deltaky = 2 * pi / params_oper.Ly
deltakx = 2 * pi / params_oper.Lx
order = spectra["dims_order"]
KX = deltakx * spectra[f"K{order[-1]}_adim"]
if self.nb_dim == 3:
deltaky = 2 * pi / params_oper.Ly
deltakz = 2 * pi / params_oper.Lz
KY = deltaky * spectra[f"K{order[1]}_adim"]
KH = np.sqrt(KX ** 2 + KY ** 2)
deltakh = max(deltakx, deltaky)
khmax_spectra = min(KX.max(), KY.max())
del KY
else:
# in 2d, vertical (here "z") is y
deltakz = 2 * pi / params_oper.Ly
KH = abs(KX)
deltakh = deltakx
khmax_spectra = KX.max()
KZ = deltakz * spectra[f"K{order[0]}_adim"]
KY = deltaky * spectra[f"K{order[1]}_adim"]
KX = deltakx * spectra[f"K{order[2]}_adim"]
KH = np.sqrt(KX ** 2 + KY ** 2)
kz_spectra = np.arange(0, KZ.max() + 1e-15, deltakz)
deltakh = max(deltakx, deltaky)
khmax_spectra = min(KX.max(), KY.max())
nkh_spectra = max(2, int(khmax_spectra / deltakh))
kh_spectra = deltakh * np.arange(nkh_spectra)
del KY
# get one-sided frequencies
omegas = spectra["omegas"]
nomegas = omegas.size // 2 + 1
......@@ -727,11 +735,17 @@ class SpatioTemporalSpectraNS:
del spectra
# total kinetic energy
spectra_kzkhomega["spectrum_K"] = 0.5 * (
spectra_kzkhomega["spectrum_vx"]
+ spectra_kzkhomega["spectrum_vy"]
+ spectra_kzkhomega["spectrum_vz"]
)
if self.nb_dim == 3:
spectra_kzkhomega["spectrum_K"] = 0.5 * (
spectra_kzkhomega["spectrum_vx"]
+ spectra_kzkhomega["spectrum_vy"]
+ spectra_kzkhomega["spectrum_vz"]
)
else:
spectra_kzkhomega["spectrum_K"] = 0.5 * (
spectra_kzkhomega["spectrum_ux"]
+ spectra_kzkhomega["spectrum_uy"]
)
# potential energy
try:
......
......@@ -28,49 +28,46 @@ A1 = "float[:]"
# @boost
def compute_spectrum_kzkhomega(
field_k0k1k2omega: A3, khs: A1, kzs: A1, KX: A3f64, KZ: A3f64, KH: A3f64
field_k0k1omega: A3, khs: A1, kzs: A1, KX: A3f64, KZ: A3f64, KH: A3f64
):
"""Compute the kz-kh-omega spectrum."""
deltakh = khs[1]
deltakz = kzs[1]
raise NotImplementedError
nkh = len(khs)
nkz = len(kzs)
nk0, nk1, nk2, nomega = field_k0k1k2omega.shape
nk0, nk1, nomega = field_k0k1omega.shape
spectrum_kzkhomega = np.zeros(
(nkz, nkh, nomega), dtype=field_k0k1k2omega.dtype
(nkz, nkh, nomega), dtype=field_k0k1omega.dtype
)
for ik0 in range(nk0):
for ik1 in range(nk1):
for ik2 in range(nk2):
values = field_k0k1k2omega[ik0, ik1, ik2, :]
kx = KX[ik0, ik1, ik2]
if kx != 0.0:
# warning: we should also consider another condition
# (kx != kx_max) but it is not necessary here mainly
# because of dealiasing
values = 2 * values
kappa = KH[ik0, ik1, ik2]
ikh = int(kappa / deltakh)
kz = abs(KZ[ik0, ik1, ik2])
ikz = int(round(kz / deltakz))
if ikz >= nkz - 1:
ikz = nkz - 1
if ikh >= nkh - 1:
ikh = nkh - 1
for i, value in enumerate(values):
spectrum_kzkhomega[ikz, ikh, i] += value
else:
coef_share = (kappa - khs[ikh]) / deltakh
for i, value in enumerate(values):
spectrum_kzkhomega[ikz, ikh, i] += (
1 - coef_share
) * value
spectrum_kzkhomega[ikz, ikh + 1, i] += coef_share * value
values = field_k0k1omega[ik0, ik1, :]
kx = KX[ik0, ik1]
if kx != 0.0:
# warning: we should also consider another condition
# (kx != kx_max) but it is not necessary here mainly
# because of dealiasing
values = 2 * values
kappa = KH[ik0, ik1]
ikh = int(kappa / deltakh)
kz = abs(KZ[ik0, ik1])
ikz = int(round(kz / deltakz))
if ikz >= nkz - 1:
ikz = nkz - 1
if ikh >= nkh - 1:
ikh = nkh - 1
for i, value in enumerate(values):
spectrum_kzkhomega[ikz, ikh, i] += value
else:
coef_share = (kappa - khs[ikh]) / deltakh
for i, value in enumerate(values):
spectrum_kzkhomega[ikz, ikh, i] += (
1 - coef_share
) * value
spectrum_kzkhomega[ikz, ikh + 1, i] += coef_share * value
# get one-sided spectrum in the omega dimension
nomega = nomega // 2 + 1
......@@ -84,7 +81,7 @@ def compute_spectrum_kzkhomega(
def _sum_wavenumber2D(field, KX, kx_max):
n0, n1 = field.shape[:-1]
n0, n1 = field.shape[:2]
result = 0.0
for i0 in range(n0):
for i1 in range(n1):
......
......@@ -242,21 +242,36 @@ class TestForcingOutput(TestSimulBase):
energy_K_mean = means["EK"][:-1].mean()
sim2.output.temporal_spectra.load_time_series()
tspectra_mean = sim2.output.temporal_spectra.compute_spectra()
omegas = tspectra_mean["omegas"]
spectra_omega = sim2.output.temporal_spectra.compute_spectra()
omegas = spectra_omega["omegas"]
delta_omega = omegas[1]
sim2.output.temporal_spectra.plot_spectra()
spatiotemporal_spectra = sim2.output.spatiotemporal_spectra
series_kxky = spatiotemporal_spectra.load_time_series()
spectra_kxkykzomega = spatiotemporal_spectra.compute_spectra()
spectra_kxkyomega = spatiotemporal_spectra.compute_spectra()
spectra_omega_from_spatiotemp = (
spatiotemporal_spectra.compute_temporal_spectra()
)
deltakx = 2 * np.pi / self.params.oper.Lx
order = spectra_kxkyomega["dims_order"]
KX = deltakx * spectra_kxkyomega[f"K{order[1]}_adim"]
kx_max = self.params.oper.nx // 2 * deltakx
assert kx_max == KX.max()
from fluidsim.solvers.ns2d.output.spatiotemporal_spectra import (
_sum_wavenumber2D,
)
def sum_wavenumber(field):
return _sum_wavenumber2D(field, KX, kx_max)
tspectrum_mean = (
tspectra_mean["spectrum_ux"] + tspectra_mean["spectrum_uy"]
spectra_omega["spectrum_ux"] + spectra_omega["spectrum_uy"]
)
energy_tspect_mean = 0.5 * delta_omega * tspectrum_mean.sum()
......@@ -265,6 +280,64 @@ class TestForcingOutput(TestSimulBase):
energy_K_mean / energy_tspect_mean
)
spectra_kzkhomega = spatiotemporal_spectra.save_spectra_kzkhomega()
delta_kz = spectra_kzkhomega["kz_spectra"][1]
delta_kh = spectra_kzkhomega["kh_spectra"][1]
delta_omega = spectra_kzkhomega["omegas"][1]
coef = delta_kz * delta_kh * delta_omega
for letter in "xy":
vi_fft = series_kxky[f"u{letter}_Fourier"]
spectrum_kxkykzomega = spectra_kxkyomega["spectrum_u" + letter]
spectrum_omega = spectra_omega["spectrum_u" + letter]
spectrum_omega_from_spatiotemp = spectra_omega_from_spatiotemp[
"spectrum_u" + letter
]
spectrum_kzkhomega = spectra_kzkhomega["spectrum_u" + letter]
E_series_kxkykz = 0.5 * sum_wavenumber(
(abs(vi_fft) ** 2).mean(axis=-1)
)
assert E_series_kxkykz > 0, (letter, vi_fft)
E_kxkykzomega = (
0.5
* delta_omega
* sum_wavenumber(spectrum_kxkykzomega.sum(axis=-1))
)
E_omega = 0.5 * delta_omega * spectrum_omega.sum()
E_omega_from_spatiotemp = (
0.5 * delta_omega * spectrum_omega_from_spatiotemp.sum()
)
E_kzkhomega = 0.5 * coef * spectrum_kzkhomega.sum()
assert np.allclose(E_omega, E_kxkykzomega), (
letter,
E_kxkykzomega / E_omega,
)
assert np.allclose(E_series_kxkykz, E_kxkykzomega), (
letter,
E_kxkykzomega / E_series_kxkykz,
)
assert np.allclose(E_series_kxkykz, E_kzkhomega), (
letter,
E_kzkhomega / E_series_kxkykz,
)
assert np.allclose(E_omega, E_omega_from_spatiotemp), (
letter,
E_omega,
E_omega_from_spatiotemp,
)
sim2.output.spatiotemporal_spectra.plot_kzkhomega(
key_field="ux", equation="kh=1"
)
sim2.output.increments.load()
sim2.output.increments.plot()
sim2.output.increments.load_pdf_from_file()
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment