diff --git a/doc/examples/internal_waves_coriolis/simul_idempotent.py b/doc/examples/internal_waves_coriolis/simul_idempotent.py index 8d48202c8..1ebe9dd8a 100644 --- a/doc/examples/internal_waves_coriolis/simul_idempotent.py +++ b/doc/examples/internal_waves_coriolis/simul_idempotent.py @@ -112,7 +112,7 @@ where $C$ is a constant of order 1. """ - n = 2 + n = 4 C = 1.0 dx = lx / nx U = amplitude * omega_f diff --git a/doc/examples/simul_ns3d_forced_isotropic.py b/doc/examples/simul_ns3d_forced_isotropic.py index 8ff4b4cf1..da541c2f0 100644 --- a/doc/examples/simul_ns3d_forced_isotropic.py +++ b/doc/examples/simul_ns3d_forced_isotropic.py @@ -14,7 +14,7 @@ help="Number of grid points in the x direction.", ) parser.add_argument( - "--t_end", type=float, default=8.0, help="End time of the simulation" + "--t_end", type=float, default=80.0, help="End time of the simulation" ) parser.add_argument( "--order", @@ -34,7 +34,7 @@ params = Simul.create_default_params() -params.output.sub_directory = "examples" +params.output.sub_directory = "examples/test_iso3d" ny = nz = nx Lx = 3 @@ -77,8 +77,8 @@ params.output.periods_save.spatial_means = 0.1 params.output.periods_save.spectra = 0.1 params.output.periods_save.spect_energy_budg = 0.1 +params.output.periods_save.kolmo_law = 0.1 params.output.spectra.kzkh_periodicity = 1 - sim = Simul(params) sim.time_stepping.start() diff --git a/fluidsim/base/output/base.py b/fluidsim/base/output/base.py index 4fedb8ef4..9867e463c 100644 --- a/fluidsim/base/output/base.py +++ b/fluidsim/base/output/base.py @@ -617,7 +617,7 @@ def _create_file_from_dict_arrays( if os.path.exists(path_file): print("file NOT created since it already exists!") elif mpi.rank == 0: - + with h5py.File(path_file, "w") as file: file.attrs["date saving"] = str(datetime.datetime.now()).encode() file.attrs["name_solver"] = self.output.name_solver diff --git a/fluidsim/base/output/kolmo_law3d.py b/fluidsim/base/output/kolmo_law3d.py index 32dad68d5..f3e38a3dc 100644 --- a/fluidsim/base/output/kolmo_law3d.py +++ b/fluidsim/base/output/kolmo_law3d.py @@ -147,7 +147,7 @@ def __init__(self, output): "drz": rz_store, "dr": r_store, } - self.pow_store = 5/4 + self.pow_store = 5 / 4 pow_store = self.pow_store for i in range(n_store): index = ((i + 1) / n_store) ** (pow_store) @@ -411,7 +411,7 @@ def plot_kolmo_law( ax2.set_xscale("log") if slope is not None: ax2.plot(posxprime, check_slope, label="slope") - # ax2.plot(posx, E_k / (rad ** (coef_comp2)), label="4*E") + # ax2.plot(posx, E_k / (rad ** (coef_comp2)), label="4*E") ax2.plot(posx, pos2ycomp, label=f"S2 compense {coef_comp2:.2g}") ax2.set_title(f"tmin={tmin:.2g}, tmax={tmax:.2g}") @@ -470,15 +470,15 @@ def compute(self): J_A_xyz = [None] * 3 E_k_mean = 0.0 - K_k = np.ones_like(K) + K_k = np.ones_like(K) E_k_proc = np.mean(K) collect_E_k = mpi.comm.gather(E_k_proc, root=0) if mpi.rank == 0: E_k_mean = np.mean(collect_E_k) else: - E_k_mean = None + E_k_mean = None - E_k_mean = mpi.comm.bcast(E_k_mean, root = 0) + E_k_mean = mpi.comm.bcast(E_k_mean, root=0) E_k = E_k_mean * K_k val = None @@ -499,7 +499,7 @@ def compute(self): for ind_j in range(3): tmp += 4 * tf_vi[ind_j] * tf_vjvi[ind_i, ind_j].conj() - tmp = 1j * tmp.imag + tmp = 1j * tmp.imag J_K_xyz[ind_i] = self.sim.oper.ifft(tmp) if "b" in keys_state_phys: @@ -507,8 +507,7 @@ def compute(self): S2_K_xyz = 2 * E_k - 2 * self.sim.oper.ifft(val) -# S2_K_xyz = 2 * K - S2_K_xyz - + # S2_K_xyz = 2 * K - S2_K_xyz J_k_z = np.zeros([n_store, n_store]) J_k_h = np.zeros([n_store, n_store]) @@ -564,16 +563,18 @@ def compute(self): if "b" in keys_state_phys: J_a_xyz = np.array(J_A_xyz) - + pow_store = self.pow_store for index, value in np.ndenumerate(J_k_xyz[2]): if "b" in keys_state_phys: rh_i = floor( - ((self.rhrz["rh"][index] / self.rh_max) ** (1/pow_store)) * n_store + ((self.rhrz["rh"][index] / self.rh_max) ** (1 / pow_store)) + * n_store ) rz_i = floor( - ((self.rhrz["rz"][index] / self.rh_max) ** (1/pow_store)) * n_store + ((self.rhrz["rz"][index] / self.rh_max) ** (1 / pow_store)) + * n_store ) count[rh_i, rz_i] += 1 J_a_hv_prov["J_a_z"][rh_i, rz_i] += J_a_xyz[2][index] @@ -583,7 +584,8 @@ def compute(self): else: r_i = floor( - ((self.rhrz["r"][index] / self.r_max) ** (1/pow_store)) * n_store + ((self.rhrz["r"][index] / self.r_max) ** (1 / pow_store)) + * n_store ) count_iso[r_i] += 1 J_k_xyz_prov["J_k_prov"][r_i] += J_k_xyz_pro[index] diff --git a/fluidsim/base/output/spectra.py b/fluidsim/base/output/spectra.py index ab1bc9447..41ca33202 100644 --- a/fluidsim/base/output/spectra.py +++ b/fluidsim/base/output/spectra.py @@ -165,6 +165,7 @@ def _plot_ndim( coef_compensate=0, coef_plot_k3=None, coef_plot_k53=None, + coef_plot_k43=None, coef_plot_k2=None, xlim=None, ylim=None, @@ -261,6 +262,10 @@ def _plot_ndim( to_plot = coef_plot_k2 * ks_no0 ** (-2) * coef_norm ax.plot(ks, to_plot, "k:", label=r"$\propto k^{-2}$") + if coef_plot_k43 is not None: + to_plot = coef_plot_k43 * ks_no0 ** (-4.0 / 3) * coef_norm + ax.plot(ks, to_plot, "k:", label=r"$\propto k^{-4/3}$") + if xlim is not None: ax.set_xlim(xlim) diff --git a/fluidsim/solvers/ns3d/output/spectra.py b/fluidsim/solvers/ns3d/output/spectra.py index 6b864a9ab..c003f81ce 100644 --- a/fluidsim/solvers/ns3d/output/spectra.py +++ b/fluidsim/solvers/ns3d/output/spectra.py @@ -267,8 +267,8 @@ def _plot_times( to_plot = coef_plot_k53 * ks_no0 ** (-5.0 / 3) * coef_norm ax.plot(ks[1:], to_plot[1:], "k-.") - if coef_plot_k2 is not None: - to_plot = coef_plot_k2 * ks_no0 ** (-2) * coef_norm + if coef_plot_k1 is not None: + to_plot = coef_plot_k1 * ks_no0 ** (-2) * coef_norm ax.plot(ks[1:], to_plot[1:], "k--") if xlim is not None: @@ -287,6 +287,7 @@ def plot1d( coef_plot_k3=None, coef_plot_k53=None, coef_plot_k2=None, + coef_plot_k43=None, xlim=None, ylim=None, directions=None, @@ -302,6 +303,7 @@ def plot1d( coef_plot_k3=coef_plot_k3, coef_plot_k53=coef_plot_k53, coef_plot_k2=coef_plot_k2, + coef_plot_k43=coef_plot_k43, xlim=xlim, ylim=ylim, ndim=1, diff --git a/fluidsim/solvers/ns3d/test_solver.py b/fluidsim/solvers/ns3d/test_solver.py index 96a95c25b..dc978f3d3 100644 --- a/fluidsim/solvers/ns3d/test_solver.py +++ b/fluidsim/solvers/ns3d/test_solver.py @@ -120,7 +120,7 @@ def init_params(cls): periods[key] = 0.1 params.output.spectra.kzkh_periodicity = 2 - + Lx, Ly, Lz = params.oper.Lx, params.oper.Ly, params.oper.Lz nx, ny, nz = params.oper.nx, params.oper.ny, params.oper.nz probes_region = (0.0, Lx, 0.0, Ly, 0.0, Lz) @@ -633,4 +633,4 @@ def test_milestone(self): if __name__ == "__main__": unittest.main() -#mpirun -n 2 python -m pytest fluidsim/solvers/ns3d/test_solver.py::TestOutput::test_output -v -s +# mpirun -n 2 python -m pytest fluidsim/solvers/ns3d/test_solver.py::TestOutput::test_output -v -s