spatiotemporal_spectra.py 9.47 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
"""
Spatiotemporal Spectra (:mod:`fluidsim.solvers.ns3d.output.spatiotemporal_spectra`)
===================================================================================

Provides:

.. autoclass:: SpatioTemporalSpectraNS3D
   :members:
   :private-members:

"""

Jason Reneuve's avatar
Jason Reneuve committed
13
14
15
16
17
from pathlib import Path

import numpy as np
import h5py

18
19
20
21
from fluidsim.base.output.spatiotemporal_spectra import SpatioTemporalSpectra


class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
Jason Reneuve's avatar
Jason Reneuve committed
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
    def loop_spectra_kzkhomega(self, spectrum_k0k1k2omega, khs, KH, kzs, KZ):
        """Compute the kz-kh-omega spectrum."""
        deltakh = khs[1]
        deltakz = kzs[1]
        nkh = len(khs)
        nkz = len(kzs)
        nk0, nk1, nk2, nomega = spectrum_k0k1k2omega.shape
        spectrum_kzkhomega = np.zeros((nkz, nkh, nomega))
        for ik0 in range(nk0):
            for ik1 in range(nk1):
                for ik2 in range(nk2):
                    value = spectrum_k0k1k2omega[ik0, ik1, ik2, :]
                    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
                        spectrum_kzkhomega[ikz, ikh, :] += value
                    else:
                        coef_share = (kappa - khs[ikh]) / deltakh
                        spectrum_kzkhomega[ikz, ikh, :] += (
                            1 - coef_share
                        ) * value
                        spectrum_kzkhomega[ikz, ikh + 1, :] += coef_share * value
Jason Reneuve's avatar
Jason Reneuve committed
49
50
51
52
53
54
55
56
57
58

        # get one-sided spectrum in the omega dimension
        nomega = nomega // 2 + 1
        spectrum_onesided = np.zeros((nkz, nkh, nomega))
        spectrum_onesided[:, :, 0] = spectrum_kzkhomega[:, :, 0]
        spectrum_onesided[:, :, 1:] = 0.5 * (
            spectrum_kzkhomega[:, :, 1:nomega]
            + spectrum_kzkhomega[:, :, -1:-nomega:-1]
        )
        return spectrum_onesided / (deltakz * deltakh)
Jason Reneuve's avatar
Jason Reneuve committed
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

    def save_spectra_kzkhomega(self, tmin=0, tmax=None):
        """save the spatiotemporal spectra, with a cylindrical average in k-space"""
        if tmax is None:
            tmax = self.sim.params.time_stepping.t_end

        # compute spectra
        dict_spectra = self.compute_spectra(tmin=tmin, tmax=tmax)

        # get kz, kh
        oper = self.sim.oper
        order = dict_spectra["dims_order"]
        KZ = oper.deltakz * dict_spectra[f"K{order[0]}_adim"]
        KY = oper.deltaky * dict_spectra[f"K{order[1]}_adim"]
        KX = oper.deltaky * dict_spectra[f"K{order[1]}_adim"]
        KH = np.sqrt(KX ** 2 + KY ** 2)

        kz_spectra = np.arange(0, KZ.max() + 1e-15, oper.deltakz)

        deltakh = oper.deltakh
        khmax_spectra = min(KX.max(), KY.max())
        nkh_spectra = max(2, int(khmax_spectra / deltakh))
        kh_spectra = deltakh * np.arange(nkh_spectra)

Jason Reneuve's avatar
Jason Reneuve committed
83
84
85
86
87
        # get one-sided frequencies
        omegas = dict_spectra["omegas"]
        nomegas = omegas.size // 2 + 1
        omegas_onesided = abs(omegas[:nomegas])

Jason Reneuve's avatar
Jason Reneuve committed
88
89
90
91
        # perform cylindrical average
        dict_spectra_kzkhomega = {
            "kz_spectra": kz_spectra,
            "kh_spectra": kh_spectra,
Jason Reneuve's avatar
Jason Reneuve committed
92
            "omegas": omegas_onesided,
Jason Reneuve's avatar
Jason Reneuve committed
93
94
95
96
97
98
99
100
        }
        for key, data in dict_spectra.items():
            if not key.startswith("spect"):
                continue
            dict_spectra_kzkhomega[key] = self.loop_spectra_kzkhomega(
                data, kh_spectra, KH, kz_spectra, KZ
            )

101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
        # total kinetic/potential energy
        dict_spectra_kzkhomega["spect_K"] = 0.5 * (
            dict_spectra_kzkhomega["spect_vx"]
            + dict_spectra_kzkhomega["spect_vy"]
            + dict_spectra_kzkhomega["spect_vz"]
        )

        # potential energy
        try:
            N = self.sim.params.N
            dict_spectra_kzkhomega["spect_A"] = (
                0.5 / N ** 2 * dict_spectra_kzkhomega["spect_b"]
            )
        except AttributeError:
            pass

Jason Reneuve's avatar
Jason Reneuve committed
117
118
119
        # save to file
        path_file = Path(self.sim.output.path_run) / "spatiotemporal_spectra.h5"
        with h5py.File(path_file, "w") as file:
120
121
            file.attrs["tmin"] = tmin
            file.attrs["tmax"] = tmax
Jason Reneuve's avatar
Jason Reneuve committed
122
123
124
            for key, val in dict_spectra_kzkhomega.items():
                file.create_dataset(key, data=val)

Jason Reneuve's avatar
Jason Reneuve committed
125
126
127
128
129
130
131
132
133
134
    def plot_kzkhomega(
        self,
        key_field=None,
        tmin=0,
        tmax=None,
        equation=None,
        cmap=None,
        vmin=None,
        vmax=None,
    ):
Jason Reneuve's avatar
Jason Reneuve committed
135
        """plot the spatiotemporal spectra, with a cylindrical average in k-space"""
Jason Reneuve's avatar
Jason Reneuve committed
136
137
        if key_field is None:
            key_field = self.keys_fields[0]
Jason Reneuve's avatar
Jason Reneuve committed
138
139
        if tmax is None:
            tmax = self.sim.params.time_stepping.t_end
Jason Reneuve's avatar
Jason Reneuve committed
140
141
        if cmap is None:
            cmap = "viridis"
Jason Reneuve's avatar
Jason Reneuve committed
142

Jason Reneuve's avatar
Jason Reneuve committed
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
        path_file = path_file = (
            Path(self.sim.output.path_run) / "spatiotemporal_spectra.h5"
        )

        dict_spectra_kzkhomega = {}

        key_spect = "spect_" + key_field

        # load spectra from file if it exists
        if path_file.is_file():
            # we should check if times match?
            print("loading spectra from file...")
            with h5py.File(path_file, "r") as file:
                for key in file.keys():
                    dict_spectra_kzkhomega[key] = file[key][...]
        else:
            # compute spectra and save to file, then load
            self.save_spectra_kzkhomega(tmin=tmin, tmax=tmax)
            with h5py.File(path_file, "r") as file:
                for key in file.keys():
                    dict_spectra_kzkhomega[key] = file[key][...]

        # slice along equation
        if equation is None:
            equation = f"omega=0"
        if equation.startswith("omega="):
            omega = eval(equation[len("omega=") :])
            omegas = dict_spectra_kzkhomega["omegas"]
            iomega = abs(omegas - omega).argmin()
            spect = dict_spectra_kzkhomega[key_spect][:, :, iomega]
            xaxis = dict_spectra_kzkhomega["kh_spectra"]
            yaxis = dict_spectra_kzkhomega["kz_spectra"]
            xlabel = r"$k_h$"
            ylabel = r"$k_z$"
            omega = omegas[iomega]
            equation = r"$\omega=$" + f"{omega:.2g}"
179
180
181
182
183
184
            # use reduced frequency for stratified fluids
            try:
                N = self.sim.params.N
                equation = r"$\omega/N=$" + f"{omega/N:.2g}"
            except AttributeError:
                pass
Jason Reneuve's avatar
Jason Reneuve committed
185
186
187
188
189
        elif equation.startswith("kh="):
            kh = eval(equation[len("kh=") :])
            kh_spectra = dict_spectra_kzkhomega["kh_spectra"]
            ikh = abs(kh_spectra - kh).argmin()
            spect = dict_spectra_kzkhomega[key_spect][:, ikh, :].transpose()
190

Jason Reneuve's avatar
Jason Reneuve committed
191
192
            xaxis = dict_spectra_kzkhomega["kz_spectra"]
            yaxis = dict_spectra_kzkhomega["omegas"]
193
194
195
196
197
198
199
            # use reduced frequency for stratified fluids
            try:
                N = self.sim.params.N
                yaxis /= N
            except AttributeError:
                pass

Jason Reneuve's avatar
Jason Reneuve committed
200
            xlabel = r"$k_z$"
201
            ylabel = r"$\omega/N$"
Jason Reneuve's avatar
Jason Reneuve committed
202
203
204
205
206
207
208
            kh = kh_spectra[ikh]
            equation = r"$k_h=$" + f"{kh:.2g}"
        elif equation.startswith("kz="):
            kz = eval(equation[len("kz=") :])
            kz_spectra = dict_spectra_kzkhomega["kz_spectra"]
            ikz = abs(kz_spectra - kz).argmin()
            spect = dict_spectra_kzkhomega[key_spect][ikz, :, :].transpose()
209

Jason Reneuve's avatar
Jason Reneuve committed
210
211
            xaxis = dict_spectra_kzkhomega["kh_spectra"]
            yaxis = dict_spectra_kzkhomega["omegas"]
212
213
214
215
216
217
218
            # use reduced frequency for stratified fluids
            try:
                N = self.sim.params.N
                yaxis /= N
            except AttributeError:
                pass

Jason Reneuve's avatar
Jason Reneuve committed
219
            xlabel = r"$k_h$"
220
            ylabel = r"$\omega/N$"
Jason Reneuve's avatar
Jason Reneuve committed
221
222
223
224
225
226
227
228
229
230
231
232
            kz = kz_spectra[ikz]
            equation = r"$k_z=$" + f"{kz:.2g}"
        else:
            raise NotImplementedError(
                "equation must start with 'omega=', 'kh=' or 'kz='"
            )

        # plot
        fig, ax = self.output.figure_axe()
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)

Jason Reneuve's avatar
Jason Reneuve committed
233
234
235
236
237
238
239
240
        if vmin is None:
            vmin = np.log10(spect[np.isfinite(spect)].min())
        if vmax is None:
            vmax = np.log10(spect[np.isfinite(spect)].max())

        im = ax.pcolormesh(
            xaxis, yaxis, np.log10(spect), cmap=cmap, vmin=vmin, vmax=vmax
        )
Jason Reneuve's avatar
Jason Reneuve committed
241
242
243
244
245
246
247
248
249
250
251
252
253
254
        fig.colorbar(im)

        ax.set_title(
            f"spatiotemporal spectra {equation}\n"
            f"tmin={tmin:.2g}, tmax={tmax:.2g}\n" + self.output.summary_simul
        )

        # add dispersion relation : omega = N * kh / sqrt(kh ** 2 + kz ** 2)
        try:
            N = self.sim.params.N
        except AttributeError:
            return
        if equation.startswith(r"$\omega"):
            kz_disp = (N ** 2 / omega ** 2 - 1) * xaxis
255
            ax.step(xaxis, kz_disp, "k", linewidth=2)
Jason Reneuve's avatar
Jason Reneuve committed
256
        elif equation.startswith(r"$k_h"):
257
            omega_disp = kh / np.sqrt(kh ** 2 + xaxis ** 2)
258
            ax.step(xaxis, omega_disp, "k", linewidth=2)
Jason Reneuve's avatar
Jason Reneuve committed
259
        elif equation.startswith(r"$k_z"):
260
            omega_disp = xaxis / np.sqrt(xaxis ** 2 + kz ** 2)
261
            ax.step(xaxis, omega_disp, "k", linewidth=2)
Jason Reneuve's avatar
Jason Reneuve committed
262
263
264
265
266
267
        else:
            raise ValueError("wrong equation for dispersion relation")

        # reset axis limits after plotting dispersion relation
        ax.set_xlim((xaxis.min(), xaxis.max()))
        ax.set_ylim((yaxis.min(), yaxis.max()))