Source code for fluidsim.solvers.plate2d.output.spatial_means

"""
Spatial means (:mod:`fluidsim.solvers.plate2d.output.spatial_means`)
==========================================================================


Provides:

.. autoclass:: SpatialMeansPlate2D
   :members:
   :private-members:

"""

import os
import numpy as np


from fluiddyn.util import mpi

from fluidsim.base.output.spatial_means import SpatialMeansBase


[docs] class SpatialMeansPlate2D(SpatialMeansBase): r"""Compute, save, load and plot spatial means. .. |p| mathmacro:: \partial If only :math:`W` is forced and dissipated, the energy budget is quite simple and can be written as: .. math:: \p_t E_w = - C_{w\rightarrow z} - C_{w\rightarrow \chi} + P_w - D_w, \p_t E_z = + C_{w\rightarrow z}, \p_t E_\chi = + C_{w\rightarrow \chi}, where .. math:: C_{w\rightarrow z} = \sum_{\mathbf{k}} k^4\mathcal{R}(\hat w \hat z^*), C_{w\rightarrow \chi} = -\sum_{\mathbf{k}} \mathcal{R}( \widehat{\{ w, z\}} \hat \chi ^* ), P_w = \sum_{\mathbf{k}} \mathcal{R}( \hat f_w \hat w^* ) and .. math:: D_w = 2 \nu_\alpha \sum_{\mathbf{k}} k^{2\alpha} E_w(k). """ def _save_one_time(self): tsim = self.sim.time_stepping.t self.t_last_save = tsim ( Ek_fft, El_fft, Ee_fft, conversion_k_to_l_fft, conversion_l_to_e_fft, ) = self.output.compute_energies_conversion_fft() energy_k = self.sum_wavenumbers(Ek_fft) energy_l = self.sum_wavenumbers(El_fft) energy_e = self.sum_wavenumbers(Ee_fft) conversion_k_to_l = self.sum_wavenumbers(conversion_k_to_l_fft) conversion_l_to_e = self.sum_wavenumbers(conversion_l_to_e_fft) f_d, f_d_hypo = self.sim.compute_freq_diss() epsK = self.sum_wavenumbers(f_d[0] * 2 * Ek_fft) epsK_hypo = self.sum_wavenumbers(f_d_hypo[0] * 2 * Ek_fft) # assert that z is not dissipated assert not (f_d[1].any() and f_d_hypo[1].any()) if self.sim.params.forcing.enable: deltat = self.sim.time_stepping.deltat state_spect = self.sim.state.state_spect w_fft = state_spect.get_var("w_fft") forcing_fft = self.sim.forcing.get_forcing() Fw_fft = forcing_fft.get_var("w_fft") # assert that z in not forced Fz_fft = forcing_fft.get_var("z_fft") assert np.allclose( abs(Fz_fft).max(), 0.0 ), f"abs(Fz_fft).max(): {abs(Fz_fft).max()}" P1_fft = np.real(w_fft.conj() * Fw_fft) P2_fft = (abs(Fw_fft) ** 2) * deltat / 2 P1 = self.sum_wavenumbers(P1_fft) P2 = self.sum_wavenumbers(P2_fft) if mpi.rank == 0: energy = energy_k + energy_l + energy_e epsK_tot = epsK + epsK_hypo self.file.write(f"####\ntime = {tsim:17.13f}\n") self.file.write( f"E = {energy:31.26e} ; E_k = {energy_k:11.6e} ; " f"E_l = {energy_l:11.6e} ; E_e = {energy_e:11.6e}\n" f"epsK = {epsK:11.6e} ; epsK_hypo = {epsK_hypo:11.6e} ; " f"epsK_tot = {epsK + epsK_hypo:11.6e}\n" f"conv_k_to_l = {conversion_k_to_l:11.6e} : conv_l_to_e = {conversion_l_to_e:11.6e}\n" ) if self.sim.params.forcing.enable: self.file.write( f"P1 = {P1:11.6e} ; P2 = {P2:11.6e} ; P_tot = {P1 + P2:11.6e} \n" ) self.file.flush() os.fsync(self.file.fileno()) if self.has_to_plot and mpi.rank == 0: self.ax_a.plot(tsim, energy, "k.") self.ax_a.plot(tsim, energy_k, "r.") self.ax_a.plot(tsim, energy_l, "b.") self.ax_a.plot(tsim, energy_e, "y.") self.axe_b.plot(tsim, epsK_tot, "k.") self.axe_b.plot(tsim, conversion_k_to_l, "c.") self.axe_b.plot(tsim, conversion_l_to_e, "g.") if self.sim.params.forcing.enable: self.axe_b.plot(tsim, P1 + P2, "m.") if tsim - self.t_last_show >= self.period_show: self.t_last_show = tsim fig = self.ax_a.get_figure() fig.canvas.draw() def load(self): dict_results = {"name_solver": self.output.name_solver} with open(self.path_file) as file_means: lines = file_means.readlines() lines_t = [] lines_E = [] lines_P = [] lines_epsK = [] for il, line in enumerate(lines): if line.startswith("time ="): lines_t.append(line) if line.startswith("E ="): lines_E.append(line) if line.startswith("P1 ="): lines_P.append(line) if line.startswith("epsK ="): lines_epsK.append(line) nt = len(lines_t) t = np.empty(nt) E = np.empty(nt) E_k = np.empty(nt) E_l = np.empty(nt) E_e = np.empty(nt) if self.sim.params.forcing.enable: P1 = np.empty(nt) P2 = np.empty(nt) P_tot = np.empty(nt) epsK = np.empty(nt) epsK_hypo = np.empty(nt) epsK_tot = np.empty(nt) for il in range(nt): line = lines_t[il] words = line.split() t[il] = float(words[2]) line = lines_E[il] words = line.split() E[il] = float(words[2]) E_k[il] = float(words[6]) E_l[il] = float(words[10]) E_e[il] = float(words[14]) if self.sim.params.forcing.enable: line = lines_P[il] words = line.split() P1[il] = float(words[2]) P2[il] = float(words[6]) P_tot[il] = float(words[10]) line = lines_epsK[il] words = line.split() epsK[il] = float(words[2]) epsK_hypo[il] = float(words[6]) epsK_tot[il] = float(words[10]) dict_results["t"] = t dict_results["E"] = E dict_results["E_k"] = E_k dict_results["E_l"] = E_l dict_results["E_e"] = E_e if self.sim.params.forcing.enable: dict_results["P1"] = P1 dict_results["P2"] = P2 dict_results["P_tot"] = P_tot dict_results["epsK"] = epsK dict_results["epsK_hypo"] = epsK_hypo dict_results["epsK_tot"] = epsK_tot return dict_results def plot(self, with_dtE=False): dict_results = self.load() t = dict_results["t"] E = dict_results["E"] E_k = dict_results["E_k"] E_l = dict_results["E_l"] E_e = dict_results["E_e"] epsK = dict_results["epsK"] epsK_hypo = dict_results["epsK_hypo"] epsK_tot = dict_results["epsK_tot"] if with_dtE: nt = len(t) dtE = np.empty(nt) for il in range(nt - 2): dtE[il + 1] = (E[il + 2] - E[il]) / (t[il + 2] - t[il]) dtE[0] = (E[1] - E[0]) / (t[1] - t[0]) dtE[nt - 1] = (E[nt - 1] - E[nt - 2]) / (t[nt - 1] - t[nt - 2]) width_axe = 0.85 height_axe = 0.8 x_left_axe = 0.12 z_bottom_axe = 0.15 size_axe = [x_left_axe, z_bottom_axe, width_axe, height_axe] fig, ax1 = self.output.figure_axe(size_axe=size_axe) ax1.set_xlabel("t") ax1.set_ylabel("Energies") ax1.plot(t, E, "k", linewidth=2) ax1.plot(t, E_k, "r", linewidth=2) ax1.plot(t, E_l, "b", linewidth=2) ax1.plot(t, E_e, "y", linewidth=2) size_axe[1] = z_bottom_axe fig, ax1 = self.output.figure_axe(size_axe=size_axe) ax1.set_xlabel("$t$") ax1.set_ylabel(r"$\varepsilon_{tot}$, $P_{tot}$, $\partial_t E$") ax1.plot(t, epsK, "k--", linewidth=2) ax1.plot(t, epsK_hypo, "k:", linewidth=2) ax1.plot(t, epsK_tot, "k", linewidth=2) if self.sim.params.forcing.enable: P_tot = dict_results["P_tot"] ax1.plot(t, P_tot, "m", linewidth=2) if with_dtE: ax1.plot(t, dtE, "b", linewidth=2) should_be_zeros = dtE + epsK_tot if self.sim.params.forcing.enable: should_be_zeros -= P_tot ax1.plot(t, should_be_zeros, "g", linewidth=2)