Commit 5497c80e authored by Pierre Augier's avatar Pierre Augier
Browse files

Fix test E_omega_from_spatiotemp

parent 65939f474ae6
Pipeline #20730 failed with stage
in 2 minutes and 5 seconds
......@@ -147,6 +147,13 @@ class TemporalSpectra3D(SpecificOutput):
self.file_max_size = params_tspec.file_max_size
self.SAVE_AS_FLOAT32 = params_tspec.SAVE_AS_FLOAT32
if self.SAVE_AS_FLOAT32:
size_1_number = 4e-6
self.datatype = np.float32
else:
size_1_number = 8e-6
self.datatype = np.float64
if self.nb_dim == 3:
X, Y, Z = oper.get_XYZ_loc()
else:
......@@ -287,11 +294,6 @@ class TemporalSpectra3D(SpecificOutput):
if self.probes_nb_loc > 0:
self._init_new_file(tmin_file=self.sim.time_stepping.t)
if self.SAVE_AS_FLOAT32:
size_1_number = 4e-6
else:
size_1_number = 8e-6
# size of a single write: nb_fields * probes_nb_loc + time
probes_write_size = (
len(self.keys_fields) * self.probes_nb_loc + 1
......@@ -339,6 +341,7 @@ class TemporalSpectra3D(SpecificOutput):
f"probes_{key}_loc",
(self.probes_nb_loc, 1),
maxshape=(self.probes_nb_loc, None),
dtype=self.datatype,
)
def _write_to_file(self, data):
......
......@@ -80,6 +80,20 @@ def compute_spectrum_kzkhomega(
return spectrum_onesided / (deltakz * deltakh)
def _sum_wavenumber3D(field, KX, kx_max):
n0, n1, n2 = field.shape[:3]
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
class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
def _get_path_saved_spectra(self, tmin, tmax, dtype, save_urud):
base = f"periodogram_{tmin}_{tmax}"
......@@ -488,26 +502,15 @@ class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
order = spectra["dims_order"]
KX = spectra[f"K{order[2]}_adim"]
kx_max = self.sim.params.oper.nx // 2 * 2 * pi / self.sim.params.oper.Lx
def _sum_wavenumber(field):
n0, n1, n2 = field.shape[:3]
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
deltakx = 2 * pi / self.sim.params.oper.Lx
kx_max = self.sim.params.oper.nx // 2 * deltakx
print(f"in compute_temporal_spectra {kx_max = }")
# average over Fourier space (kx,ky,kz)
for key, spectrum in spectra.items():
if not key.startswith("spectrum_"):
continue
tspectra[key] = _sum_wavenumber(spectrum)
tspectra[key] = _sum_wavenumber3D(spectrum, KX, kx_max)
tspectra["omegas"] = spectra["omegas"]
......
......@@ -118,7 +118,7 @@ class TestOutput(TestSimulBase):
params.output.temporal_spectra.probes_deltax = Lx / nx
params.output.temporal_spectra.probes_deltay = Ly / ny
params.output.temporal_spectra.probes_deltaz = Lz / nz
params.output.temporal_spectra.SAVE_AS_FLOAT32 = True
params.output.temporal_spectra.SAVE_AS_FLOAT32 = False
nx, ny, nz = params.oper.nx, params.oper.ny, params.oper.nz
params.output.spatiotemporal_spectra.probes_region = (
......@@ -126,7 +126,7 @@ class TestOutput(TestSimulBase):
ny // 2,
nz // 2,
)
params.output.spatiotemporal_spectra.SAVE_AS_COMPLEX64 = True
params.output.spatiotemporal_spectra.SAVE_AS_COMPLEX64 = False
@pytest.mark.filterwarnings("ignore:divide by zero encountered in log10")
def test_output(self):
......@@ -259,18 +259,14 @@ class TestOutput(TestSimulBase):
KX = deltakx * spectra_kxkykzomega[f"K{order[2]}_adim"]
kx_max = self.params.oper.nx // 2 * deltakx
assert kx_max == KX.max()
from fluidsim.solvers.ns3d.output.spatiotemporal_spectra import (
_sum_wavenumber3D,
)
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
return _sum_wavenumber3D(field, KX, kx_max)
spectra_kzkhomega = spatiotemporal_spectra.save_spectra_kzkhomega(
save_urud=True
......@@ -301,6 +297,10 @@ class TestOutput(TestSimulBase):
* 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()
# `:-1` because the last time is saved twice in spatial_means
E_mean = means["E" + letter][:-1].mean()
......@@ -325,6 +325,15 @@ class TestOutput(TestSimulBase):
E_series_kxkykz / E_mean,
)
assert np.allclose(E_omega, E_omega_from_spatiotemp), (
letter,
E_omega,
E_omega_from_spatiotemp,
)
# print(f"{spectrum_omega.sum() / spectrum_omega_from_spatiotemp.sum() = }")
# this condition is not exactly fulfilled (why?)
# assert np.allclose(
# spectrum_omega, spectrum_omega_from_spatiotemp
# ), (
......@@ -332,16 +341,9 @@ class TestOutput(TestSimulBase):
# spectrum_omega,
# spectrum_omega_from_spatiotemp,
# spectrum_omega / spectrum_omega_from_spatiotemp,
# spectrum_omega.sum() / spectrum_omega_from_spatiotemp.sum(),
# )
assert (
spectrum_omega.sum() == spectrum_omega_from_spatiotemp.sum()
), (
letter,
spectrum_omega.sum(),
spectrum_omega_from_spatiotemp.sum(),
)
spectrum_Khd = spectra_kzkhomega["spectrum_Khd"]
spectrum_vz = spectra_kzkhomega["spectrum_vz"]
......
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