Read about our upcoming Code of Conduct on this issue

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

spatiotemporal_spectra: fix test energy computations

parent c906d97bb7ad
......@@ -38,7 +38,7 @@ shortlog:
@hg log -M -r$(RELEASE): --template '- {desc|firstline} (:rev:`{node|short}`)\n'
black:
black -l 82 fluidsim scripts bench doc lib
black -l 82 fluidsim scripts bench doc lib --exclude __pythran__
tests:
pytest -v lib
......
......@@ -550,10 +550,11 @@ class SpatioTemporalSpectra(SpecificOutput):
fs=self.f_sample,
scaling="spectrum",
detrend=False,
return_onesided=False,
)
return freq, spectrum / self.domega
def compute_spectra(self, tmin=0, tmax=None, dtype=None):
def compute_temporal_spectra(self, tmin=0, tmax=None, dtype=None):
"""compute spatiotemporal spectra from files"""
if tmax is None:
tmax = self.sim.params.time_stepping.t_end
......
......@@ -26,8 +26,8 @@ A1 = "float[:]"
@boost
def loop_spectra_kzkhomega(
spectrum_k0k1k2omega: A4, khs: A1, KH: A3, kzs: A1, KZ: A3
def compute_spectrum_kzkhomega(
field_k0k1k2omega: A4, khs: A1, kzs: A1, KX: A3, KZ: A3, KH: A3
):
"""Compute the kz-kh-omega spectrum."""
deltakh = khs[1]
......@@ -35,15 +35,20 @@ def loop_spectra_kzkhomega(
nkh = len(khs)
nkz = len(kzs)
nk0, nk1, nk2, nomega = spectrum_k0k1k2omega.shape
nk0, nk1, nk2, nomega = field_k0k1k2omega.shape
spectrum_kzkhomega = np.zeros(
(nkz, nkh, nomega), dtype=spectrum_k0k1k2omega.dtype
(nkz, nkh, nomega), dtype=field_k0k1k2omega.dtype
)
for ik0 in range(nk0):
for ik1 in range(nk1):
for ik2 in range(nk2):
values = spectrum_k0k1k2omega[ik0, ik1, ik2, :]
values = field_k0k1k2omega[ik0, ik1, ik2, :]
kx = KX[ik0, ik1, ik2]
# TODO: other condition (kx != kx_max)
if kx != 0.0:
values = 2 * values
kappa = KH[ik0, ik1, ik2]
ikh = int(kappa / deltakh)
kz = abs(KZ[ik0, ik1, ik2])
......@@ -66,7 +71,7 @@ def loop_spectra_kzkhomega(
nomega = nomega // 2 + 1
spectrum_onesided = np.zeros((nkz, nkh, nomega))
spectrum_onesided[:, :, 0] = spectrum_kzkhomega[:, :, 0]
spectrum_onesided[:, :, 1:] = 0.5 * (
spectrum_onesided[:, :, 1:] = (
spectrum_kzkhomega[:, :, 1:nomega]
+ spectrum_kzkhomega[:, :, -1:-nomega:-1]
)
......@@ -91,7 +96,7 @@ class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
# compute spectra
print("Computing spectra...")
spectra = self.compute_spectra(tmin=tmin, tmax=tmax, dtype=dtype)
spectra = self.compute_temporal_spectra(tmin=tmin, tmax=tmax, dtype=dtype)
# get kz, kh
params_oper = self.sim.params.oper
......@@ -111,7 +116,7 @@ class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
nkh_spectra = max(2, int(khmax_spectra / deltakh))
kh_spectra = deltakh * np.arange(nkh_spectra)
del KX, KY
del KY
# get one-sided frequencies
omegas = spectra["omegas"]
......@@ -127,8 +132,8 @@ class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
for key, data in spectra.items():
if not key.startswith("spectrum_"):
continue
spectra_kzkhomega[key] = loop_spectra_kzkhomega(
np.ascontiguousarray(data), kh_spectra, KH, kz_spectra, KZ
spectra_kzkhomega[key] = compute_spectrum_kzkhomega(
np.ascontiguousarray(data), kh_spectra, kz_spectra, KX, KZ, KH
)
del spectra
......@@ -166,8 +171,8 @@ class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
for key, data in spectra.items():
if not key.startswith("spectrum_"):
continue
spectra_urud_kzkhomega[key] = loop_spectra_kzkhomega(
np.ascontiguousarray(data), kh_spectra, KH, kz_spectra, KZ
spectra_urud_kzkhomega[key] = compute_spectrum_kzkhomega(
np.ascontiguousarray(data), kh_spectra, kz_spectra, KX, KZ, KH
)
spectra_kzkhomega[key] = spectra_urud_kzkhomega[key]
......@@ -348,10 +353,6 @@ class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
keys=("vx", "vy"), tmin=tmin, tmax=tmax, dtype=dtype
)
# get the sampling frequency
times = series["times"]
f_sample = 1 / np.mean(times[1:] - times[:-1])
# toroidal/poloidal decomposition
# urx_fft, ury_fft contain shear modes!
vx_fft = series["vx_Fourier"]
......
......@@ -43,10 +43,10 @@ class TestSimulBase(TestSimul):
Lx = 6.0
params.oper.Lx = Lx
params.oper.Ly = Lx * 3 // 4
params.oper.Ly = Lx * params.oper.ny / params.oper.nx
try:
params.oper.nz = nx // 2
params.oper.Lz = Lx // 2
params.oper.Lz = Lx * params.oper.nz / params.oper.nx
except AttributeError:
pass
......@@ -239,6 +239,26 @@ class TestOutput(TestSimulBase):
spatiotemporal_spectra = sim3.output.spatiotemporal_spectra
series = spatiotemporal_spectra.load_time_series()
tspectra = spatiotemporal_spectra.compute_temporal_spectra()
deltakx = 2 * pi / self.params.oper.Lx
order = tspectra["dims_order"]
KX = deltakx * tspectra[f"K{order[2]}_adim"]
kx_max = self.params.oper.nx // 2 * deltakx
def sum_wavenumber(field):
n0, n1, n2 = field.shape
result = 0.0
for i0 in range(n0):
for i1 in range(n1):
for i2 in range(n2):
value = field[i0, i1, i2]
kx = KX[i0, i1, i2]
if kx != 0.0 and kx != kx_max:
value *= 2
result += value
return result
spectra = spatiotemporal_spectra.save_spectra_kzkhomega(
save_urud=True
)
......@@ -248,21 +268,29 @@ class TestOutput(TestSimulBase):
delta_omega = spectra["omegas"][1]
coef = delta_kz * delta_kh * delta_omega
print(spectra["kz_spectra"])
print(spectra["kh_spectra"])
for letter in "xyz":
vi_fft = series[f"v{letter}_Fourier"]
tspectrum = tspectra["spectrum_v" + letter]
spectrum_vi = spectra["spectrum_v" + letter]
# TODO: compute energy from vi_fft and spectrum_vi
energy_fft = (0.5 * abs(vi_fft) ** 2).mean(axis=-1).sum()
assert energy_fft > 0, (letter, vi_fft)
energy_serie = 0.5 * sum_wavenumber(
(abs(vi_fft) ** 2).mean(axis=-1)
)
assert energy_serie > 0, (letter, vi_fft)
energy_tspect = (
0.5 * delta_omega * sum_wavenumber(tspectrum.sum(axis=-1))
)
energy_spe = 0.5 * coef * spectrum_vi.sum()
# TODO: fix this and plug this condition
assert np.allclose(energy_fft, energy_spe), (
assert np.allclose(energy_serie, energy_tspect), (
letter,
energy_tspect / energy_serie,
)
assert np.allclose(energy_serie, energy_spe), (
letter,
energy_spe / energy_fft,
energy_spe / energy_serie,
)
spectrum_Khd = spectra["spectrum_Khd"]
......@@ -271,9 +299,8 @@ class TestOutput(TestSimulBase):
# because k \cdot \hat v = 0, for kz = 0, Khd = 0
assert np.allclose(spectrum_Khd[0].sum(), 0.0)
# DONE: understand why we need the `1:`
# energy in the mode kx=ky=kz=0 is not zero, for any field in state_spect.
# because k \cdot \hat v = 0, for kh = 0, Kz = 0
# `1:` because the energy in the mode kx=ky=kz=0 is not zero.
assert np.allclose(spectrum_vz[1:, 0, :].sum(), 0.0)
sim3.output.spatiotemporal_spectra.plot_kzkhomega(
......
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