Modif Kolmolaw aniso

Author: agoua

doc/examples/simul_ns3d_forced_isotropic.py

@@ -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()

fluidsim/base/output/base.py

@@ -617,7 +617,6 @@ 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

fluidsim/base/output/kolmo_law3d.py

@@ -66,7 +66,7 @@ class KolmoLaw(SpecificOutput):
 
     This output saves the components in the spherical basis of the vectors
     :math:`\J_\alpha` averaged over the azimutal angle (i.e. as a function of
-    :math:`r_h` and :math:`r_z`).
+    :math:`r_h` and :math:`r_v`).
 
     We can take the example of the quantity :math:`\langle | \delta b |^2
     \delta \v \rangle_\x` to explain how these quantities are computed. Using
@@ -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)
@@ -418,7 +418,7 @@ def plot_kolmo_law(
                     ax2.set_xscale("log")
                 if slope is not None:
                     ax2.plot(posxprime, check_slope, label=f"$r^{2}$")
- #               ax2.plot(posx, E_k / (rad ** (coef_comp2)), label="4*E")
+
                 ax2.plot(posx, pos2ycomp)
                 ax2.set_title(f"tmin={tmin:.2g}, tmax={tmax:.2g}")
 
@@ -496,15 +496,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
@@ -525,16 +525,18 @@ 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
+            mom = 1j * tmp.imag
 
             J_K_xyz[ind_i] = self.sim.oper.ifft(tmp)
+
             if "b" in keys_state_phys:
                 J_A_xyz[ind_i] = self.sim.oper.ifft(mom)
 
         S2_K_xyz = 2 * E_k - 2 * self.sim.oper.ifft(val)
 
-#        S2_K_xyz = 2 * K - S2_K_xyz
-
+        nh_store = n_store
+        nv_store = n_store
 
         J_k_v = np.zeros([nh_store, nv_store])
         J_k_h = np.zeros([nh_store, nv_store])
@@ -592,7 +594,7 @@ 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]):  # average on each process
@@ -622,29 +624,24 @@ def compute(self):
 
             if "b" in keys_state_phys:
                 J_a_hv_average["J_a_v"][rh_i, rv_i] += J_a_xyz[2][index]
-                    ((self.rhrz["rh"][index] / self.rh_max) ** (1/pow_store)) * n_store
-                )
-                    ((self.rhrz["rz"][index] / self.rh_max) ** (1/pow_store)) * n_store
-
-            else:
-                r_i = floor(
-                    ((self.rhrz["r"][index] / self.r_max) ** (1/pow_store)) * n_store
+                J_a_hv_average["J_a_h"][rh_i, rv_i] += np.sqrt(
+                    J_a_xyz[0] ** 2 + J_a_xyz[1] ** 2
                 )
-                count_iso[r_i] += 1
-                J_k_xyz_prov["J_k_prov"][r_i] += J_k_xyz_pro[index]
-                S2_k_xyz_prov["S2_k_prov"][r_i] += S2_k_xyz[index]
 
-        if mpi.nb_proc == 1:
-            J_k_hv_prov["count"] = J_a_hv_prov["count"] = count
+        if mpi.nb_proc > 0:  # average on one process
 
-        if mpi.nb_proc > 0:
-
-            collect_J_k_prov = mpi.comm.gather(J_k_xyz_prov["J_k_prov"], root=0)
-            collect_S2_k_prov = mpi.comm.gather(
-                S2_k_xyz_prov["S2_k_prov"], root=0
+            collect_J_k_average = mpi.comm.gather(
+                J_k_xyz_average["J_k_average"], root=0
+            )  # gather all results on one process
+            collect_S2_k_average = mpi.comm.gather(
+                S2_k_xyz_average["S2_k_average"], root=0
             )
             collect_count_iso = mpi.comm.gather(count_iso, root=0)
 
+            collect_J_k_v = mpi.comm.gather(J_k_hv_average["J_k_v"], root=0)
+            collect_J_k_h = mpi.comm.gather(J_k_hv_average["J_k_h"], root=0)
+            collect_count = mpi.comm.gather(count, root=0)
+
             if "b" in keys_state_phys:
                 collect_J_a_v = mpi.comm.gather(J_a_hv_average["J_a_v"], root=0)
                 collect_J_a_h = mpi.comm.gather(J_a_hv_average["J_a_h"], root=0)
@@ -672,20 +669,24 @@ def compute(self):
                 J_k_xyz_average["count"] = count_final_iso
                 S2_k_xyz_average["count2"] = count_final_iso
 
-                    count_final_iso = np.sum(collect_count_iso, axis=0)
-                    J_k_xyz_prov["count"] = count_final_iso
-                    S2_k_xyz_prov["count2"] = count_final_iso
-                if "b" in keys_state_phys:
-                    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]
+                for index, value in np.ndenumerate(J_k_hv_average["J_k_v"]):
+                    if count_final[index] == 0:
+                        J_k_hv_average["J_k_v"][index] = 0.0
+                        J_k_hv_average["J_k_h"][index] = 0.0
+                        if "b" in keys_state_phys:
+                            J_a_hv_average["J_a_v"][index] = 0.0
+                            J_a_hv_average["J_a_h"][index] = 0.0
+                    else:
+                        J_k_hv_average["J_k_v"][index] = (
+                            value / count_final[index]
+                        )
+                        J_k_hv_average["J_k_h"][index] = (
+                            J_k_hv_average["J_k_h"][index] / count_final[index]
+                        )
+                        if "b" in keys_state_phys:
+                            J_a_hv_average["J_a_v"][index] = (
+                                J_a_hv_average["J_a_v"][index]
+                                / count_final[index]
                             )
                             J_a_hv_average["J_a_h"][index] = (
                                 J_a_hv_average["J_a_h"][index]

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)
 

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:

fluidsim/solvers/ns3d/test_solver.py

@@ -118,7 +118,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)
@@ -628,4 +628,3 @@ 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