Modif Kolmolaw aniso

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

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

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

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,16 +499,15 @@ 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:
                 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
-
+        #        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]

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:
@@ -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,

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