From 2e988248cd985212191bc22c86efdad3428657cd Mon Sep 17 00:00:00 2001 From: agoua Date: Thu, 16 May 2024 11:12:07 +0200 Subject: [PATCH] Modification Kolmolaw --- .../jean_zay/install/4_clone_fluid.sh | 0 .../clusters/licallo/install/4_clone_fluid.sh | 0 .../licallo/install/install_p3dfft.sh | 0 .../clusters/licallo/install/install_pfft.sh | 0 fluidsim/base/output/base.py | 1 + fluidsim/base/output/kolmo_law3d.py | 292 ++++++++++++++++-- fluidsim/solvers/ns3d/test_solver.py | 10 +- 7 files changed, 270 insertions(+), 33 deletions(-) mode change 100755 => 100644 doc/examples/clusters/jean_zay/install/4_clone_fluid.sh mode change 100755 => 100644 doc/examples/clusters/licallo/install/4_clone_fluid.sh mode change 100755 => 100644 doc/examples/clusters/licallo/install/install_p3dfft.sh mode change 100755 => 100644 doc/examples/clusters/licallo/install/install_pfft.sh diff --git a/doc/examples/clusters/jean_zay/install/4_clone_fluid.sh b/doc/examples/clusters/jean_zay/install/4_clone_fluid.sh old mode 100755 new mode 100644 diff --git a/doc/examples/clusters/licallo/install/4_clone_fluid.sh b/doc/examples/clusters/licallo/install/4_clone_fluid.sh old mode 100755 new mode 100644 diff --git a/doc/examples/clusters/licallo/install/install_p3dfft.sh b/doc/examples/clusters/licallo/install/install_p3dfft.sh old mode 100755 new mode 100644 diff --git a/doc/examples/clusters/licallo/install/install_pfft.sh b/doc/examples/clusters/licallo/install/install_pfft.sh old mode 100755 new mode 100644 diff --git a/fluidsim/base/output/base.py b/fluidsim/base/output/base.py index c33901148..4fedb8ef4 100644 --- a/fluidsim/base/output/base.py +++ b/fluidsim/base/output/base.py @@ -617,6 +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 8de883443..e6063f4b5 100644 --- a/fluidsim/base/output/kolmo_law3d.py +++ b/fluidsim/base/output/kolmo_law3d.py @@ -12,55 +12,60 @@ import itertools import numpy as np - +import os from fluiddyn.util import mpi +import h5py +import matplotlib.pyplot as plt + from .base import SpecificOutput +from math import floor + """ Conversion from cartesian coordinates system to spherical coordinate system """ class OperKolmoLaw: def __init__(self, X, Y, Z, params): - self.r = np.sqrt(X**2 +Y**2 + Z**2) + self.r = np.sqrt(X**2 + Y**2 + Z**2) self.rh = np.sqrt(X**2 + Y**2) self.rz = Z -# self.sin_theta_sin_phi = Y / self.r -# self.cos_phi_sin_theta = X / self.r -# self.sin_theta = self.rh / self.r -# self.cos_theta = Z / self.r -# self.sin_phi = Y / self.rh -# self.cos_phi = X / self.rh - def average_azimutal(self, arr): + def average_azimutal(self, arr): + avg_arr = None if mpi.nb_proc == 1: - avg_arr = np.mean(arr, axis=(0,1)) + avg_arr = np.mean(arr, axis=(0, 1)) return avg_arr - local_sum = np.sum(arr, axis=(0,1)) + local_sum = np.sum(arr, axis=(0, 1)) if mpi.rank == 0: global_arr = np.zeros(self.nz) # define array to sum on all proc for rank in range(mpi.nb_proc): if mpi.rank == 0: nz_loc = self.nzs_local[rank] # define size of array on each proc + print("nz_loc " + str(nz_loc)) if rank == 0 and mpi.rank == 0: sum = local_sum # start the sum on rank 0 else: - if mpi.rank == 0: # sum made on rank 0: receive local_array of rank + # sum made on rank 0: receive local_array of rank + if ( + mpi.rank == 0 + ): sum = np.empty(nz_loc) mpi.comm.Recv(sum, source=rank, tag=42 * rank) elif mpi.rank == rank: # send the local array to 0 mpi.comm.Send(sum_local, dest=0, tag=42 * rank) if mpi.rank == 0: # construct global sum on 0 iz_start = self.izs_start[rank] + print("iz_start " + iz_start) global_array[iz_start : iz_start + nz_loc] += sum if mpi.rank == 0: - avg_arr = global_array / self.nazim + avg_arr = global_array / (params.oper.nx * params.oper.ny) return avg_arr @@ -130,7 +135,7 @@ class KolmoLaw(SpecificOutput): """ _tag = "kolmo_law" - _name_file = "kolmo_law.h5" + _name_file = _tag + ".h5" @classmethod def _complete_params_with_default(cls, params): @@ -145,31 +150,172 @@ def __init__(self, output): period_save_kolmo_law = params.output.periods_save.kolmo_law except AttributeError: period_save_kolmo_law = 0.0 + period_save_kolmo_law = 0.1 if period_save_kolmo_law != 0.0: X, Y, Z = output.sim.oper.get_XYZ_loc() self.oper_kolmo_law = OperKolmoLaw(X, Y, Z, params) + + self.rhrz = { + "rh": self.oper_kolmo_law.rh, + "rz": self.oper_kolmo_law.rz + } + self.rh_max = np.sqrt(params.oper.Lx**2 + params.oper.Ly**2) + self.rz_max = params.oper.Lz + # n_sort = params.n_sort + self.n_sort = 50 + n_sort = self.n_sort + rh_sort = np.empty([n_sort]) + rz_sort = np.empty([n_sort]) + self.drhrz = { + "drh": rh_sort, + "drz": rz_sort, + } + + + for i in range(n_sort): + self.drhrz["drh"][i] = self.rh_max * (i + 1) / n_sort + self.drhrz["drz"][i] = self.rz_max * (i + 1) / n_sort arrays_1st_time = { - "rh": self.oper_kolmo_law.rh, - "rz": self.oper_kolmo_law.rz, + "rh_sort": self.drhrz["drh"], + "rz_sort": self.drhrz["drz"], } + else: arrays_1st_time = None + self.rhrz_sort = arrays_1st_time super().__init__( output, - period_save=period_save_kolmo_law, - arrays_1st_time=arrays_1st_time, + #period_save=period_save_kolmo_law, + period_save = params.output.periods_save.spectra, + arrays_1st_time = arrays_1st_time, ) + def _init_path_files(self): + + path_run = self.output.path_run + self.path_kolmo_law = path_run + "/kolmo_law.h5" + self.path_file = self.path_kolmo_law + + + + def _init_files(self, arrays_1st_time=None): + dict_J_k_hv = self.compute()[0] + dict_J_a_hv = self.compute()[1] + dict_J = {} + dict_J.update(dict_J_k_hv) + dict_J.update(dict_J_a_hv) + if mpi.rank == 0: +# print("dict_J_k_hv= " + str(dict_J_k_hv)) + if not os.path.exists(self.path_kolmo_law): + + self._create_file_from_dict_arrays( + self.path_kolmo_law, dict_J, arrays_1st_time + ) + self.nb_saved_times = 1 + else: + print("Fichier modif") + with h5py.File(self.path_kolmo_law, "r") as file: + dset_times = file["times"] + self.nb_saved_times = dset_times.shape[0] + 1 + print(self.nb_saved_times) + self._add_dict_arrays_to_file(self.path_kolmo_law, dict_J) + + self.t_last_save = self.sim.time_stepping.t + + def _online_save(self): + """Save the values at one time.""" + tsim = self.sim.time_stepping.t + if tsim - self.t_last_save >= self.period_save: + self.t_last_save = tsim + dict_J_k_hv = self.compute()[0] + dict_J_a_hv = self.compute()[1] + dict_J = {} + dict_J.update(dict_J_k_hv) + dict_J.update(dict_J_a_hv) + if mpi.rank == 0: + self._add_dict_arrays_to_file(self.path_kolmo_law, dict_J) + self.nb_saved_times += 1 + + def load(self): + path_file = self.path_kolmo_law + J_hv = { + "J_k_h" : None, + "J_k_z" : None, + "J_a_h" : None, + "J_a_z" : None, + "times" : None, + "count" : None + } + file = h5py.File(path_file, "r") + for key in file.keys(): + J_hv[key] = file[key] + return J_hv + + def plot_kolmo_law(self): + result = self.load() + times=result["times"] + J_k_h = result["J_k_h"][0] + J_k_z = result["J_k_z"][0] + J_a_h = result["J_a_h"][0] + J_a_z = result["J_a_z"][0] + J_tot_h = J_a_h + J_k_h + J_tot_z = J_a_z + J_k_z + count = result["count"] + + print(str(result.items())) + print("J_k_z = " + str(J_k_z[:])) + print("J_k_h = " + str(J_k_h[:])) + print("J_a_z = " + str(J_a_z[:])) + print("J_a_h = " + str(J_a_h[:])) + print("count = " + str(count[:])) + print("count_tot = " + str(320*320*80) + " " + "count_sum= " + + str(sum(sum(count[0])))) + + + posy = result["rz_sort"][:] + posx= result["rh_sort"][:] + U,V= np.meshgrid(posx,posy) + toty = J_tot_z + totx = J_tot_h + + bx = J_a_h + by = J_a_z + + kx = J_k_h + ky = J_k_z + + if mpi.rank == 0: + fig, ax = self.output.figure_axe() + self.ax = ax + ax.set_xlabel("$rh$") + ax.set_ylabel("$rz$") + ax.set_title("J_tot") + ax.quiver(posx,posy,totx,toty) + + fig2, ax2 = self.output.figure_axe() + self.ax2 = ax2 + ax2.set_xlabel("$rh$") + ax2.set_ylabel("$rz$") + ax2.set_title("J_A") + ax2.quiver(posx,posy,bx,by) + + fig3, ax3 = self.output.figure_axe() + self.ax3 = ax3 + ax3.set_xlabel("$rh$") + ax3.set_ylabel("$rz$") + ax3.set_title("J_K") + ax3.quiver(posx,posy,kx,ky) + def compute(self): """compute the values at one time.""" state = self.sim.state + params = self.sim.params state_phys = state.state_phys state_spect = state.state_spect - keys_state_phys=state.keys_state_phys + keys_state_phys = state.keys_state_phys fft = self.sim.oper.fft - letters = "xyz" tf_vi = [state_spect.get_var(f"v{letter}_fft") for letter in letters] @@ -178,11 +324,11 @@ def compute(self): tf_K = None if "b" in keys_state_phys: - b= check_strat(state_phys.get_var("b")) - tf_b = check_strat(state_spect.get_var("b_fft")) + b = state_phys.get_var("b") + tf_b = state_spect.get_var("b_fft") b2 = b * b tf_b2 = fft(b2) - tf_bv=[None] * 3 + tf_bv = [None] * 3 bv = [item * b for item in vel] for index in range(len(bv)): tf_bv[index] = fft(bv[index]) @@ -204,31 +350,117 @@ def compute(self): J_K_xyz = [None] * 3 J_A_xyz = [None] * 3 + for ind_i in range(3): tmp = tf_vi[ind_i] * tf_K.conj() if "b" in keys_state_phys: - mom = 4 * tf_bv[ind_i].conj() * tf_b + 2 * tf_b2.conj() * tf_vi[ind_i] + mom = ( + 4 * tf_bv[ind_i].conj() * tf_b + + 2 * tf_b2.conj() * tf_vi[ind_i] + ) mom.real = 0.0 for ind_j in range(3): tmp += tf_vi[ind_j] * tf_vjvi[ind_i, ind_j].conj() + tmp.real = 0.0 + J_K_xyz[ind_i] = 4 * self.sim.oper.ifft(tmp) if "b" in keys_state_phys: J_A_xyz[ind_i] = self.sim.oper.ifft(mom) - result = {} - keys = ["r", "theta", "phi"] - for index, key in enumerate(keys): - result["J_K_" + key] = self.oper_kolmo_law.average_azimutal( - J_K_xyz[index] + + n_sort = self.n_sort + J_k_z = np.zeros([n_sort, n_sort]) + J_k_h = np.zeros([n_sort, n_sort]) + + + J_a_z = np.zeros([n_sort, n_sort]) + J_a_h = np.zeros([n_sort, n_sort]) + + count_final = np.empty([n_sort, n_sort]) + + J_k_hv_prov = { + "J_k_h": J_k_h, + "J_k_z": J_k_z, + "count" : count_final, + } + J_a_hv_prov = { + "J_a_h": J_a_h, + "J_a_z": J_a_z, + "count" : count_final, + } + count = np.zeros([n_sort, n_sort], dtype = int) + rhrz_sort = self.rhrz_sort + J_k_xyz = np.array(J_K_xyz) + J_a_xyz = np.array(J_A_xyz) + + + + for index, value in np.ndenumerate(J_k_xyz[2]): + i = floor(self.rhrz["rh"][index] * n_sort / self.rh_max) + j = floor(self.rhrz["rz"][index] * n_sort / self.rz_max) + count[i, j] += 1 + J_k_hv_prov["J_k_z"][i, j] += np.abs(value) + J_k_hv_prov["J_k_h"][i, j] += np.sqrt( + J_k_xyz[0][index] ** 2 + J_k_xyz[1][index] ** 2 ) + J_a_hv_prov["J_a_z"][i, j] += np.abs(J_a_xyz[2][index]) + + J_a_hv_prov["J_a_h"][i, j] += np.sqrt( + J_a_xyz[0][index] ** 2 + J_a_xyz[1][index] ** 2 + ) + + + if mpi.nb_proc == 1: + J_k_hv_prov["count"] = J_a_hv_prov["count"] =count + + + if mpi.nb_proc > 0: + collect_J_k_z = mpi.comm.gather(J_k_hv_prov["J_k_z"], root=0) + collect_J_k_h = mpi.comm.gather(J_k_hv_prov["J_k_h"], root=0) + collect_J_a_z = mpi.comm.gather(J_a_hv_prov["J_a_z"], root=0) + collect_J_a_h = mpi.comm.gather(J_a_hv_prov["J_a_h"], root=0) + collect_count = mpi.comm.gather(count, root=0) + + + if mpi.rank == 0: + J_k_hv_prov["J_k_z"] = np.sum(collect_J_k_z, axis=0) + J_k_hv_prov["J_k_h"] = np.sum(collect_J_k_h, axis=0) + + J_a_hv_prov["J_a_z"] = np.sum(collect_J_a_z, axis=0) + J_a_hv_prov["J_a_h"] = np.sum(collect_J_a_h, axis=0) + + count_final = np.sum(collect_count, axis = 0) + J_k_hv_prov["count"] = J_a_hv_prov["count"] = count_final + + for index, value in np.ndenumerate(J_k_hv_prov["J_k_z"]): + if count_final[index] == 0: + J_k_hv_prov["J_k_z"][index] = 0.0 + J_k_hv_prov["J_k_h"][index] = 0.0 + J_a_hv_prov["J_a_z"][index] = 0.0 + J_a_hv_prov["J_a_h"][index] = 0.0 + else: + J_k_hv_prov["J_k_z"][index] = value / count_final[index] + J_k_hv_prov["J_k_h"][index] = ( + J_k_hv_prov["J_k_h"][index] / count_final[index] + ) + J_a_hv_prov["J_a_z"][index] = ( + J_a_hv_prov["J_a_z"][index] / count_final[index] + ) + J_a_hv_prov["J_a_h"][index] = ( + J_a_hv_prov["J_a_h"][index] / count_final[index] + ) + + result = J_k_hv_prov, J_a_hv_prov return result + def check_diff_methods(self): first_method = compute(self) second_method = compute_alt(self) if not np.allclose(first_method, second_method): raise RuntimeError( - "Both methods are inconsistent: " " ({self.sim.time_stepping.it = })" + "Both methods are inconsistent: " + " ({self.sim.time_stepping.it = })" ) diff --git a/fluidsim/solvers/ns3d/test_solver.py b/fluidsim/solvers/ns3d/test_solver.py index 2e480d612..96a95c25b 100644 --- a/fluidsim/solvers/ns3d/test_solver.py +++ b/fluidsim/solvers/ns3d/test_solver.py @@ -47,7 +47,7 @@ def _init_grid(params, nx): @classmethod def init_params(cls): params = cls.params = cls.Simul.create_default_params() - + print("Coucou") params.short_name_type_run = "test" params.output.sub_directory = "unittests" cls._init_grid(params, nx=cls.nx) @@ -63,7 +63,7 @@ def init_params(cls): params.oper.coef_dealiasing = 2.0 / 3 params.nu_4 = 2.0 params.nu_8 = 2.0 - + params.n_sort = 5.0 params.time_stepping.t_end = 1.5 * params.time_stepping.deltat_max params.init_fields.type = "noise" @@ -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) @@ -407,6 +407,8 @@ def customize(result, sim): assert df.loc[0, "foo"] > 0 sim.output.get_mean_values(customize=customize) + sim.output.kolmo_law.compute() + class TestInitInScript(TestSimulBase): @classmethod @@ -630,3 +632,5 @@ 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