Skip to content
Snippets Groups Projects
Commit 3e8ba6bb authored by calpe's avatar calpe
Browse files

Working on k-omega spectra.

parent 7e606cdf
No related branches found
No related tags found
No related merge requests found
......@@ -80,6 +80,12 @@
}
classes.spect_energy_budg._set_attribs(attribs)
attribs = {
"module_name": base_name_mod + ".spatio_temporal_spectra",
"class_name": "SpatioTempSpectra",
}
classes._set_child("spatio_temporal_spectra", attribs=attribs)
def compute_energies_fft(self):
"""Compute the kinetic and potential energy (k)"""
rot_fft = self.sim.state.state_spect.get_var("rot_fft")
......
......@@ -22,6 +22,7 @@
import h5py
import matplotlib.pyplot as plt
from math import pi
from fluiddyn.util import mpi
from fluiddyn.calcul.easypyfft import FFTW1DReal2Complex, FFTW2DReal2Complex
......@@ -74,9 +75,6 @@
n1 = len(
list(range(0, output.sim.oper.shapeX_loc[1], self.coef_decimate))
)
print("n0", n0)
print("n1", n1)
# 3D array (time, x, y) and init FFTW object
# self.spatio_temp = np.empty([self.nb_times_compute, n0, n1])
self.spatio_temp = np.empty([self.nb_times_compute, n0, n1 // 2 + 1])
# Initialization operators FFT and hanning windows.
self.oper_fft2 = FFTW2DReal2Complex(n0, n1)
......@@ -82,5 +80,5 @@
self.oper_fft2 = FFTW2DReal2Complex(n0, n1)
self.oper_fft1 = FFTW1DReal2Complex(self.spatio_temp.shape, axis=0)
self.nb_omegas = self.oper_fft1.shapeK[0]
self.oper_fft1 = FFTW1DReal2Complex(
(self.nb_times_compute, n0, n1 // 2 + 1), axis=0)
self.hamming = np.hanning(self.nb_times_compute)
......@@ -85,8 +83,17 @@
self.hamming = np.hanning(self.nb_times_compute)
# Compute kxs and kys with the decimate values
deltakx = 2 * np.pi / self.sim.oper.Lx
self.kxs_decimate = np.arange(0, deltakx * (n0 / 2) + deltakx, deltakx)
# 3D array
self.spatio_temp = np.empty(
[self.nb_times_compute, n0, n1 // 2 + 1], dtype=complex)
# Array omegas
deltat = self.sim.time_stepping.deltat
nt = self.nb_times_compute
print("deltat", deltat)
time_tot = deltat * nt
self.delta_omega = 2 * pi / time_tot
self.omegas = np.arange(
0, self.delta_omega * (nt // 2 + 1), self.delta_omega)
self.nb_times_in_spatio_temp = 0
......@@ -90,42 +97,8 @@
self.nb_times_in_spatio_temp = 0
if os.path.exists(self.path_file):
with h5py.File(self.path_file, "r") as f:
link_spatio_temp_spectra = f["spatio_temp_spectra"]
self.spatio_temp_spectra = link_spatio_temp_spectra[-1]
self.periods_fill = f["periods_fill"][...]
if self.sim.time_stepping.deltat != f["deltat"][...]:
raise ValueError()
else:
self.periods_fill = params.output.periods_save.spatio_temporal_spectra
if self.periods_fill > 0:
# self.periods_fill = self.periods_fill - 1
dt_output = self.periods_fill * output.sim.time_stepping.deltat
print("dt_output = ", dt_output)
duration = self.nb_times_compute * dt_output
self.delta_omega = 2 * np.pi / duration
print("duration = ", duration)
print("delta_omega = ", self.delta_omega)
self.omegas = self.delta_omega * np.arange(self.nb_omegas)
self.omega_Nyquist = np.pi / dt_output
self.omega_dealiasing = (
self.params.oper.coef_dealiasing
* np.pi
* self.params.oper.nx
/ self.params.oper.Lx
) ** 2
if self.omega_dealiasing > self.omega_Nyquist:
print("Warning: omega_dealiasing > omega_Nyquist")
def _init_files(self, dict_arrays_1time=None):
# we can not do anything when this function is called.
pass
def _init_files2(self, spatio_temp_spectra):
......@@ -127,12 +100,8 @@
def _init_files(self, dict_arrays_1time=None):
# we can not do anything when this function is called.
pass
def _init_files2(self, spatio_temp_spectra):
time_tot = (
self.sim.time_stepping.deltat
* self.nb_times_compute
* self.periods_fill
)
omegas = 2 * np.pi / time_tot * np.arange(self.nb_omegas)
""" Initialize a file to save. """
dict_arrays_1time = {
......@@ -138,3 +107,3 @@
dict_arrays_1time = {
"omegas": omegas,
"omegas": self.omegas,
"deltat": self.sim.time_stepping.deltat,
......@@ -140,7 +109,6 @@
"deltat": self.sim.time_stepping.deltat,
"nb_times_compute": self.nb_times_compute,
"periods_fill": self.periods_fill,
}
"nb_times_compute": self.nb_times_compute}
self._create_file_from_dict_arrays(
self.path_file, spatio_temp_spectra, dict_arrays_1time
)
......@@ -152,8 +120,4 @@
itsim = int(self.sim.time_stepping.t / self.sim.time_stepping.deltat)
periods_save = self.sim.params.output.periods_save.spatio_temporal_spectra
# print('it_sim = ', itsim)
# print('it_sim - it_last_run = ', itsim - self.it_last_run)
# print('periods_save = ', periods_save - 1)
if itsim - self.it_last_run >= periods_save - 1:
......@@ -159,5 +123,4 @@
if itsim - self.it_last_run >= periods_save - 1:
# print('#save_period')
self.it_last_run = itsim
field = self.sim.state.state_phys.get_var(self.key_quantity)
field_decimate = field[::self.coef_decimate, ::self.coef_decimate]
......@@ -166,7 +129,6 @@
self.nb_times_in_spatio_temp += 1
if self.nb_times_in_spatio_temp == self.nb_times_compute:
# print('#####save_spatio_temporal...')
self.nb_times_in_spatio_temp = 0
self.t_last_save = self.sim.time_stepping.t
......@@ -170,6 +132,5 @@
self.nb_times_in_spatio_temp = 0
self.t_last_save = self.sim.time_stepping.t
# It is not the best way to apply hanning.
for i, value in enumerate(self.hamming):
self.spatio_temp[i, :, :] = value * self.spatio_temp[i, :, :]
......@@ -174,7 +135,6 @@
for i, value in enumerate(self.hamming):
self.spatio_temp[i, :, :] = value * self.spatio_temp[i, :, :]
# self.spatio_fft = self.oper_fft1.fft(
# self.hamming * self.spatio_temp)
self.spatio_fft = self.oper_fft1.fft(self.spatio_temp)
if mpi.rank == 0:
......@@ -199,37 +159,8 @@
"""
return (1 / 2.) * np.abs(spatio_temporal_fft) ** 2
def compute_spatio_temporal_spectra(self, spatio_temporal_fft):
"""Compute the spectra (kx, omega) and (ky, omega)"""
energy_fft = self._compute_energy_from_spatio_temporal_fft(
spatio_temporal_fft
)
# Axes of the spatio temporal fft.
omega_axis = 0
ky_axis = 1
kx_axis = 2
delta_kx = self.oper.deltakx
delta_ky = self.oper.deltaky
# Compute energy spectra (kx, omega).
# We multiply by 2 and 2 because there are only omega>=0 and kx>=0.
# We divide by two because the energy at the zero modes shoud be
# counted only once.
E_kx_omega = 2. * 2. * energy_fft.sum(ky_axis) / (
delta_kx * self.delta_omega
)
E_kx_omega[0, :] = E_kx_omega[0, :] / 2.
E_kx_omega[:, 0] = E_kx_omega[:, 0] / 2.
return E_kx_omega
def plot(self):
""" Plot the spatio temporal spectra. """
# Load data from file.
# Spatio_temporal_fft is an array with all spatio_temporal_fft.
# We'll need to average all of them.
def plot_frequency_spectra(self):
""" Plots the frequency spectra F(\omega). """
# Load data from file
spatio_temporal_fft = self.load()
......@@ -234,5 +165,5 @@
spatio_temporal_fft = self.load()
# Average all spatio_temporal_fft. Axis = 0
# Average over all spatio_temporal_fft 3d arrays
spatio_temporal_fft = spatio_temporal_fft.mean(axis=0)
......@@ -237,6 +168,6 @@
spatio_temporal_fft = spatio_temporal_fft.mean(axis=0)
# Compute the spatio_temporal_spectra
E_kx_omega = self.compute_spatio_temporal_spectra(spatio_temporal_fft)
# print('E_kx_omega = ', E_kx_omega)
# Compute energy from spatio_temporal_fft
energy_fft = self._compute_energy_from_spatio_temporal_fft(
spatio_temporal_fft)
......@@ -242,6 +173,4 @@
# Data grid
omegas = self.omegas
kx_grid, omegas_grid = np.meshgrid(self.kxs_decimate, omegas)
omegas = np.arange(
0, self.delta_omega * (self.nb_times_compute// 2 + 1), self.delta_omega)
......@@ -247,7 +176,3 @@
fig = plt.figure()
ax = fig.gca()
ax.set_xlabel("kx")
ax.set_ylabel("omega")
ax.set_title("E_kx_omega")
E_omega = (1 / self.delta_omega) * energy_fft.sum(1).sum(1)
......@@ -253,10 +178,9 @@
# ax.pcolor(
# kx_grid, omegas_grid, E_kx_omega, \
# vmin=E_kx_omega.min(), vmax=E_kx_omega.max())
ax.pcolor(
kx_grid, omegas_grid, E_kx_omega, vmin=1e-5, vmax=E_kx_omega.max()
)
# Plot
fig, ax = plt.subplots()
ax.set_xlabel("$\omega$")
ax.set_ylabel(r"$F(\omega)$")
ax.set_title(r"$E(\omega)$")
ax.loglog(omegas, E_omega)
plt.show()
......@@ -261,2 +185,35 @@
plt.show()
def plot_omega_kx(self):
""" Plots the frequency spectra F(\omega). """
# Load data from file
spatio_temporal_fft = self.load()
# Average over all spatio_temporal_fft 3d arrays
spatio_temporal_fft = spatio_temporal_fft.mean(axis=0)
# Compute energy from spatio_temporal_fft
energy_fft = self._compute_energy_from_spatio_temporal_fft(
spatio_temporal_fft)
delta_kx = 2 * pi / self.oper.Lx
nx_decimate = len(
list(range(0, self.sim.oper.shapeX_loc[1], self.coef_decimate)))
kxs = np.arange(0, delta_kx * (nx_decimate // 2 + 1), delta_kx)
omegas = np.arange(
0, self.delta_omega * (self.nb_times_compute// 2 + 1), self.delta_omega)
E_omega_kx = (1 / self.delta_omega) * (1 / delta_kx) * energy_fft.sum(1)
# Plot
kx_grid, omega_grid = np.meshgrid(kxs, omegas)
fig, ax = plt.subplots()
ax.set_xlabel("$\omega$")
ax.set_ylabel("$k_x$")
ax.set_title(r"$E(\omega, k_x)$")
ax.pcolor(omega_grid, kx_grid, E_omega_kx)
plt.show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment