diff --git a/fluidsim/base/output/kolmo_law3d.py b/fluidsim/base/output/kolmo_law3d.py index 2b7981a5..3dcf9387 100644 --- a/fluidsim/base/output/kolmo_law3d.py +++ b/fluidsim/base/output/kolmo_law3d.py @@ -28,7 +28,9 @@ class OperKolmoLaw: def __init__(self, X, Y, Z, params): self.r = np.sqrt(X**2 + Y**2 + Z**2) + # self.r = X self.rh = np.sqrt(X**2 + Y**2) + # self.rh = Y self.rz = np.abs(Z) self.X = X self.Y = Y @@ -53,15 +55,15 @@ class KolmoLaw(SpecificOutput): .. math:: - \bnabla \cdot \left( \J_K + \J_p \right) = -4 \left( \epsK + \epsA \right), + \bnabla \cdot \left( \Jk + \Jp \right) = -4 \left( \epsK + \epsA \right), where .. math:: - \J_K(\r) \equiv + \Jk(\r) \equiv \left\langle | \delta \v |^2 \delta \v \right\rangle_\x, \\ - \J_p(\r) \equiv + \Jp(\r) \equiv \frac{1}{N^2} \left\langle | \delta b |^2 \delta \v \right\rangle_\x. This output saves the components in the spherical basis of the vectors @@ -137,7 +139,7 @@ def __init__(self, output): params.oper.Lx**2 + params.oper.Ly**2 + params.oper.Lz**2 ) - self.n_store = 120 + self.n_store = 80 n_store = self.n_store rh_store = np.empty([n_store]) rz_store = np.empty([n_store]) @@ -151,9 +153,11 @@ def __init__(self, output): pow_store = self.pow_store for i in range(n_store): index = ((i + 1) / n_store) ** (pow_store) - self.drhrz["drh"][i] = self.rh_max * index - self.drhrz["drz"][i] = self.rz_max * index - self.drhrz["dr"][i] = self.r_max * index + # self.drhrz["drh"][i] = self.rh_max * index + # self.drhrz["drz"][i] = self.rz_max * index + self.drhrz["drh"][i] = 0.5 * self.r_max * index + self.drhrz["drz"][i] = 0.5 * self.r_max * index + self.drhrz["dr"][i] = 0.5 * self.r_max * index arrays_1st_time = { "rh_store": self.drhrz["drh"], "rz_store": self.drhrz["drz"], @@ -181,14 +185,12 @@ def _init_files(self, arrays_1st_time=None): state = self.sim.state params = self.sim.params keys_state_phys = state.keys_state_phys - if "b" in keys_state_phys: - dict_J_k_r, dict_J_k_hv, dict_J_p_hv, dict_S2_k_xyz = self.compute() - dict_J = {} - dict_J.update([dict_J_k_r, dict_J_k_hv, dict_J_p_hv, dict_S2_k_xyz]) - else: - dict_J_k_r, dict_J_k_hv, dict_S2_k_xyz = self.compute() - dict_J = {} - dict_J.update([dict_J_k_r, dict_J_k_hv, dict_S2_k_xyz]) + + dict_J = {} + result = self.compute() + for name in result: + dict_J.update(name) + if mpi.rank == 0: if not os.path.exists(self.path_kolmo_law): @@ -214,18 +216,11 @@ def _online_save(self): tsim = self.sim.time_stepping.t if tsim - self.t_last_save >= self.period_save: self.t_last_save = tsim - if "b" in keys_state_phys: - dict_J_k_r, dict_J_k_hv, dict_J_p_hv, dict_S2_k_xyz = ( - self.compute() - ) - dict_J = {} - dict_J.update( - [dict_J_k_r, dict_J_k_hv, dict_J_p_hv, dict_S2_k_xyz] - ) - else: - dict_J_k_r, dict_J_k_hv, dict_S2_k_xyz = self.compute() - dict_J = {} - dict_J.update([dict_J_k_r, dict_J_k_hv, dict_S2_k_xyz]) + result = self.compute() + dict_J = {} + for name in result: + dict_J.update(name) + if mpi.rank == 0: self._add_dict_arrays_to_file(self.path_kolmo_law, dict_J) self.nb_saved_times += 1 @@ -237,21 +232,29 @@ def load(self): keys_state_phys = state.keys_state_phys if "b" in keys_state_phys: J_hv = { - "J_k_r": None, - "J_k_h": None, - "J_k_v": None, - "S2_k": None, - "J_p_h": None, - "J_p_v": None, + "Jk_average": None, + "Jk_laverage": None, + "divJk_average": None, + "divJk_hv": None, + "Jk_h": None, + "Jk_v": None, + "Jk_l": None, + "S2_k_average": None, + "Jp_h": None, + "Jp_v": None, "times": None, "count": None, } else: J_hv = { - "J_k_r": None, - "J_k_h": None, - "J_k_v": None, - "S2_k": None, + "Jk_average": None, + "Jk__laverage": None, + "divJk_average": None, + "divJk_hv": None, + "Jk_h": None, + "Jk_v": None, + "Jk_l": None, + "S2_k_average": None, "times": None, "count": None, "count2": None, @@ -261,6 +264,9 @@ def load(self): J_hv[key] = file[key] return J_hv + def div(self, A): + return 1j * (self.oper.kx) + def plot_kolmo_law( self, tmin=None, @@ -292,22 +298,15 @@ def plot_kolmo_law( if "b" in keys_state_phys: - J_k_h = result["J_k_h"][1] - J_k_v = result["J_k_v"][2] - J_p_h = result["J_p_h"][4] - J_p_v = result["J_p_v"][5] - J_tot_h = J_p_h + J_k_h - J_tot_v = J_p_v + J_k_v + Jk_h = result["Jk_h"][1] + Jk_v = result["Jk_v"][2] + Jp_h = result["Jp_h"][4] + Jp_v = result["Jp_v"][5] + J_tot_h = Jp_h + Jk_h + J_tot_v = Jp_v + Jk_v count = result["count"] print("count = " + str(count[:])) - print( - "count_tot = " - + str(320 * 320 * 80) - + " " - + "count_sum= " - + str(sum(sum(count[0]))) - ) posy = result["rz_store"][:] posx = result["rh_store"][:] @@ -315,11 +314,11 @@ def plot_kolmo_law( toty = J_tot_v totx = J_tot_h - bx = J_p_h - by = J_p_v + bx = Jp_h + by = Jp_v - kx = J_k_h - ky = J_k_v + kx = Jk_h + ky = Jk_v if mpi.rank == 0: fig, ax = self.output.figure_axe() @@ -333,59 +332,67 @@ def plot_kolmo_law( self.ax2 = ax2 ax2.set_xlabel("$rh$") ax2.set_ylabel("$rz$") - ax2.set_title("J_p") + ax2.set_title("Jp") 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_r") + ax3.set_title("Jk_average") ax3.quiver(posx, posy, kx, ky) else: - J_k_average = result["J_k_r"][imin_plot:imax_plot] - J_k_average = -J_k_average - S2_k_average = result["S2_k"][imin_plot:imax_plot] + Jk_average = result["Jk_average"][imin_plot:imax_plot] + Jk_laverage = result["Jk_laverage"][imin_plot:imax_plot] + divJk = result["divJk_average"][imin_plot:imax_plot] + + Jk_h = result["Jk_h"][imin_plot:imax_plot] + Jk_l = result["Jk_l"][imin_plot:imax_plot] + Jk_v = result["Jk_v"][imin_plot:imax_plot] + Jk_average = -Jk_average + Jk_laverage = -Jk_laverage + S2_k_average = result["S2_k_average"][imin_plot:imax_plot] count = result["count"] count2 = result["count2"] - J_k_average = np.mean(J_k_average, axis=0) + Jk_average = np.mean(Jk_average, axis=0) + Jk_laverage = np.mean(Jk_laverage, axis=0) + divJk = np.mean(divJk, axis=0) S2_k_average = np.mean(S2_k_average, axis=0) + Jk_h = np.mean(Jk_h, axis=0) + Jk_v = np.mean(Jk_v, axis=0) + Jk_hmean = np.mean(Jk_h, axis=1) + Jk_vmean = np.mean(Jk_v, axis=0) + L = 3 n = params.oper.nx dx = L / n eta = dx - print( - "count_tot = " - + str(n * n * n) - + " " - + "count_sum= " - + str((sum(count[0]))) - + " " - + str((count[0])) - ) rad = result["r_store"][:] - print(str(rad)) + posx = rad / eta unite = posx / posx + posz = result["rz_store"][:] + rtest = posz[20] + posh = result["rh_store"][:] posxprime = rad[0:5] / eta compa = 4 / 3 * unite - pos3y = J_k_average + pos3y = Jk_average pos2y = S2_k_average pos3ycomp = pos3y / (rad ** (coef_comp3)) pos2ycomp = pos2y / (rad ** (coef_comp2)) - + pos3ybis = Jk_laverage / (rad ** (coef_comp3)) if slope is not None: check_slope = 0.01 * (posxprime**slope) E_k = 4 * 1.27 * np.ones([len(rad)]) if mpi.rank == 0: fig1, ax1 = self.output.figure_axe() self.ax1 = ax1 - ax1.set_xlabel("$r/eta$") + ax1.set_xlabel("$r/\eta$") ax1.set_ylabel("$J_L(r)$") # /(r\epsilon) if scale is None: @@ -398,12 +405,40 @@ def plot_kolmo_law( if slope is not None: ax1.plot(posxprime, check_slope, label=f"$r^{2}$") ax1.plot(posx, pos3ycomp, label="$J_{LC}(r)$") + ax1.plot(posx, pos3ybis, label="$Jbis_{LC}(r)$") ax1.set_title(f"tmin={tmin:.2g}, tmax={tmax:.2g}") # ax1.plot(posx,compa, label = "4/3") ax1.set_xscale("log") # ax1.set_yscale("log") ax1.legend() + fig5, ax5 = self.output.figure_axe() + self.ax5 = ax5 + ax5.set_xlabel("$r/eta$") + ax5.set_ylabel("$J_v(r)$") + # /(r\epsilon) + if scale is None: + ax5.set_yscale("log") + ax5.set_xscale("log") + ax5.set_ylim([0.001, 2.0]) + else: + ax5.set_yscale(f"{scale}") + ax5.set_xscale("log") + if slope is not None: + ax5.plot(posxprime, check_slope, label=f"$r^{2}$") + ax5.plot(posz / eta, -Jk_vmean / posz, label="$J_{KV}(rv)$") + ax5.plot(posh / eta, -Jk_hmean / posh, label="$J_{Kh}(rh)$") + ax5.plot( + posx, + -(Jk_hmean * posh + Jk_vmean * posz) / rad**2, + label="$J_{Kh}(rh)$", + ) + ax5.set_title(f"tmin={tmin:.2g}, tmax={tmax:.2g}") + # ax1.plot(posx,compa, label = "4/3") + ax5.set_xscale("log") + # ax1.set_yscale("log") + ax5.legend() + fig2, ax2 = self.output.figure_axe() self.ax2 = ax2 ax2.set_xlabel("$r/eta$") @@ -442,6 +477,35 @@ def plot_kolmo_law( ax3.legend() + fig4, ax4 = self.output.figure_axe() + self.ax4 = ax4 + ax4.set_xlabel("$rh$") + ax4.set_ylabel("$rz$") + ax4.set_xlim([1.2, 3.5]) + ax4.set_title(f"$Jk(r_h,r_v)$") + ax4.quiver(posh, posz, Jk_h / posh, Jk_v / posz) + ax4.set_title(f"tmin={tmin:.2g}, tmax={tmax:.2g}") + + ax4.legend() + + fig6, ax6 = self.output.figure_axe() + self.ax6 = ax6 + ax6.set_xlabel("$r/\eta$") + ax6.set_ylabel("$-div(J)/4\epsilon$") + if scale is None: + ax6.set_yscale("log") + ax6.set_xscale("log") + ax6.set_ylim([0.9, 10.0]) + else: + ax6.set_yscale(f"{scale}") + ax6.set_xscale("log") + + ax6.plot(posx, -divJk / 4) + ax6.plot(posx, np.ones([len(rad)]), label="theoretical value=1") + ax6.set_title(f"tmin={tmin:.2g}, tmax={tmax:.2g}") + + ax6.legend() + def compute(self): """compute the values at one time.""" state = self.sim.state @@ -456,13 +520,16 @@ def compute(self): letters = "xyz" n_store = self.n_store vol = params.oper.nx * params.oper.ny * params.oper.nz - + kx = self.sim.oper.Kx + ky = self.sim.oper.Ky + kz = self.sim.oper.Kz tf_vi = [state_spect.get_var(f"v{letter}_fft") for letter in letters] vel = [state_phys.get_var(f"v{letter}") for letter in letters] tf_vjvi = np.empty((3, 3), dtype=object) tf_K = None K = None + if "b" in keys_state_phys: b = state_phys.get_var("b") tf_b = state_spect.get_var("b_fft") @@ -490,24 +557,35 @@ def compute(self): vj = state_phys.get_var("v" + letter_j) tf_vjvi[ind_i, ind_j] = tf_vjvi[ind_j, ind_i] = fft(vi * vj) - J_k_r = [None] * 3 + Jk_r = [None] * 3 + Jk_r_fft = [None] * 3 if "b" in keys_state_phys: - J_p_xyz = [None] * 3 - + Jp_r = [None] * 3 + Jp_r_fft = [None] * 3 E_k_mean = 0.0 K_k = np.ones_like(K) E_k_proc = np.mean(K) + bv_proc = None + collect_bv = None + if "b" in keys_state_phys: + den_flux = np.ones_like(bv) + bv_proc = np.mean(bv, axis=1) if mpi.nb_proc > 1: collect_E_k = mpi.comm.gather(E_k_proc, root=0) + # if "b" in keys_state_phys: + # collect_bv = mpi.comm.gather(bv_proc, root=0) if mpi.rank == 0: E_k_mean = np.mean(collect_E_k) + # bv_mean = np.mean(collect_bv) else: E_k_mean = None E_k_mean = mpi.comm.bcast(E_k_mean, root=0) else: - E_k_mean = Ek_proc + E_k_mean = E_k_proc + # bv_mean = bv_proc E_k = E_k_mean * K_k + # bv_flux = bv_mean * den_flux val = None for ind_i in range(3): @@ -529,77 +607,116 @@ def compute(self): tmp = 1j * tmp.imag mom = 1j * tmp.imag - - J_k_r[ind_i] = self.sim.oper.ifft(tmp) + Jk_r_fft[ind_i] = tmp + Jk_r[ind_i] = self.sim.oper.ifft(tmp) if "b" in keys_state_phys: - J_p_xyz[ind_i] = self.sim.oper.ifft(mom) + Jp_r[ind_i] = self.sim.oper.ifft(mom) + Jp_r_fft[ind_i] = mom + S2_k_r = 2 * E_k - 2 * self.sim.oper.ifft(val) - S2_K_xyz = 2 * E_k - 2 * self.sim.oper.ifft(val) + Jk_r_fft = np.array(Jk_r_fft) + divJk_fft = 1j * (kx * Jk_r_fft[0] + ky * Jk_r_fft[1] + kz * Jk_r_fft[2]) + divJk = self.sim.oper.ifft(divJk_fft) + + if "b" in keys_state_phys: + Jp_r_fft = np.array(Jp_r_fft) + divJp_fft = 1j * ( + kx * Jp_r_fft[0] + ky * Jp_r_fft[1] + kz * Jp_r_fft[2] + ) + divJp = self.sim.oper.ifft(divJp_fft) 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]) - J_k_average = np.zeros([n_store]) + Jk_v = np.zeros([nh_store, nv_store]) + Jk_h = np.zeros([nh_store, nv_store]) + Jk_l = np.zeros([nh_store, nv_store]) + divJk_hv = np.zeros([nh_store, nv_store]) + Jk_average = np.zeros([n_store]) + divJk_average = np.zeros([n_store]) + divJp_average = np.zeros([n_store]) + Jk_laverage = np.zeros([n_store]) S2_k_average = np.zeros([n_store]) if "b" in keys_state_phys: - J_p_v = np.zeros([nh_store, nv_store]) - J_p_h = np.zeros([nh_store, nv_store]) - + Jp_v = np.zeros([nh_store, nv_store]) + Jp_h = np.zeros([nh_store, nv_store]) + KP_ex = np.zeros([nh_store, nv_store]) count_final = np.empty([nh_store, nv_store]) count_final_iso = np.empty([n_store]) - J_k_hv_average = { - "J_k_h": J_k_h, - "J_k_v": J_k_v, + Jk_hv_average = { + "Jk_h": Jk_h, + "Jk_v": Jk_v, + "Jk_l": Jk_l, + "divJk_hv": divJk_hv, "count": count_final, } - J_k_r_average = { - "J_k_r": J_k_average, + Jk_r_average = { + "Jk_average": Jk_average, + "Jk_laverage": Jk_laverage, + "divJk_average": divJk_average, "count": count_final_iso, } - S2_k_xyz_average = { - "S2_k": S2_k_average, + S2_k_r_average = { + "S2_k_average": S2_k_average, "count2": count_final_iso, } if "b" in keys_state_phys: - J_p_hv_average = { - "J_p_h": J_p_h, - "J_p_v": J_p_v, + Jp_hv_average = { + "Jp_h": Jp_h, + "Jp_v": Jp_v, + "count": count_final, + } + Jp_r_average = { + "divJp_average": divJp_average, + } + KP_exchange = { + "KP_ex": KP_ex, "count": count_final, } - count = np.zeros([nh_store, nv_store], dtype=int) count_iso = np.zeros([n_store], dtype=int) rhrz = self.rhrz_store - J_k_r = np.array(J_k_r) - S2_k_xyz = np.array(S2_K_xyz) - J_k_r_pro = np.empty_like(J_k_r) - - for index, value in np.ndenumerate( - self.rhrz["r"][:] - ): # Longitudinal projection + Jk_r = np.array(Jk_r) + S2_k_r = np.array(S2_k_r) + Jk_r_pro = np.empty_like(X) + Jk_h_pro = np.empty_like(X) + Jk_l_pro = np.empty_like(X) + + for index, value in np.ndenumerate(self.rhrz["r"][:]): + if self.rhrz["rh"][index] == 0.0: + Jk_h_pro[index] = 0.0 + else: + Jk_h_pro[index] = ( + Jk_r[0][index] * X[index] + Jk_r[1][index] * Y[index] + ) / self.rhrz["rh"][index] + # Longitudinal projection if value == 0.0: - J_k_r_pro[index] = 0.0 + Jk_r_pro[index] = 0.0 + Jk_l_pro[index] = 0.0 else: - J_k_r_pro[index] = ( - J_k_r[0][index] * X[index] - + J_k_r[1][index] * Y[index] - + J_k_r[2][index] * Z[index] + Jk_r_pro[index] = ( + Jk_r[0][index] * X[index] + + Jk_r[1][index] * Y[index] + + Jk_r[2][index] * Z[index] + ) / value + + Jk_l_pro[index] = ( + Jk_h_pro[index] * self.rhrz["rh"][index] + + Jk_r[2][index] * Z[index] ) / value if "b" in keys_state_phys: - J_p_xyz = np.array(J_p_xyz) + Jp_r = np.array(Jp_r) pow_store = self.pow_store - for index, value in np.ndenumerate(J_k_r[2]): # average on each process + for index, value in np.ndenumerate(Jk_r[2]): # average on each process rh_i = floor( ((self.rhrz["rh"][index] / self.rh_max) ** (1 / pow_store)) @@ -617,100 +734,136 @@ def compute(self): count[rh_i, rv_i] += 1 count_iso[r_i] += 1 - J_k_hv_average["J_k_v"][rh_i, rv_i] += value - J_k_hv_average["J_k_h"][rh_i, rv_i] += np.sqrt( - J_k_r[0] ** 2 + J_k_r[1] ** 2 - ) - J_k_r_average["J_k_r"][r_i] += J_k_r_pro[index] - S2_k_xyz_average["S2_k"][r_i] += S2_k_xyz[index] + Jk_hv_average["Jk_v"][rh_i, rv_i] += value + Jk_hv_average["Jk_h"][rh_i, rv_i] += Jk_h_pro[index] + Jk_hv_average["Jk_l"][rh_i, rv_i] += Jk_l_pro[index] + Jk_hv_average["divJk_hv"][rh_i, rv_i] += divJk[index] + + Jk_r_average["Jk_average"][r_i] += Jk_r_pro[index] + Jk_r_average["divJk_average"][r_i] += divJk[index] + Jk_r_average["Jk_laverage"][r_i] += Jk_l_pro[index] + S2_k_r_average["S2_k_average"][r_i] += S2_k_r[index] if "b" in keys_state_phys: - J_p_hv_average["J_p_v"][rh_i, rv_i] += J_p_xyz[2][index] - J_p_hv_average["J_p_h"][rh_i, rv_i] += np.sqrt( - J_p_xyz[0] ** 2 + J_p_xyz[1] ** 2 + Jp_hv_average["Jp_v"][rh_i, rv_i] += Jp_r[2][index] + Jp_hv_average["Jp_h"][rh_i, rv_i] += np.sqrt( + Jp_r[0][index] ** 2 + Jp_r[1][index] ** 2 ) if mpi.nb_proc > 1: # average on one process - collect_J_k_average = mpi.comm.gather( - J_k_r_average["J_k_r"], root=0 + collect_Jk_average = mpi.comm.gather( + Jk_r_average["Jk_average"], root=0 + ) + + collect_divJk_average = mpi.comm.gather( + Jk_r_average["divJk_average"], root=0 + ) + + collect_Jk_laverage = mpi.comm.gather( + Jk_r_average["Jk_laverage"], root=0 ) # gather all results on one process + collect_S2_k_average = mpi.comm.gather( - S2_k_xyz_average["S2_k"], root=0 + S2_k_r_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_Jk_v = mpi.comm.gather(Jk_hv_average["Jk_v"], root=0) + collect_Jk_h = mpi.comm.gather(Jk_hv_average["Jk_h"], root=0) + collect_Jk_l = mpi.comm.gather(Jk_hv_average["Jk_l"], root=0) + collect_divJk_hv = mpi.comm.gather(Jk_hv_average["divJk_hv"], root=0) collect_count = mpi.comm.gather(count, root=0) if "b" in keys_state_phys: - collect_J_p_v = mpi.comm.gather(J_p_hv_average["J_p_v"], root=0) - collect_J_p_h = mpi.comm.gather(J_p_hv_average["J_p_h"], root=0) + collect_Jp_v = mpi.comm.gather(Jp_hv_average["Jp_v"], root=0) + collect_Jp_h = mpi.comm.gather(Jp_hv_average["Jp_h"], root=0) if mpi.rank == 0: - J_k_r_average["J_k_r"] = np.sum(collect_J_k_average, axis=0) - S2_k_xyz_average["S2_k"] = np.sum(collect_S2_k_average, axis=0) - J_k_hv_average["J_k_v"] = np.sum(collect_J_k_v, axis=0) - J_k_hv_average["J_k_h"] = np.sum(collect_J_k_h, axis=0) + Jk_r_average["Jk_average"] = np.sum(collect_Jk_average, axis=0) + Jk_r_average["divJk_average"] = np.sum( + collect_divJk_average, axis=0 + ) + Jk_r_average["Jk_laverage"] = np.sum(collect_Jk_laverage, axis=0) + S2_k_r_average["S2_k_average"] = np.sum( + collect_S2_k_average, axis=0 + ) + Jk_hv_average["Jk_v"] = np.sum(collect_Jk_v, axis=0) + Jk_hv_average["Jk_h"] = np.sum(collect_Jk_h, axis=0) + Jk_hv_average["Jk_l"] = np.sum(collect_Jk_l, axis=0) + Jk_hv_average["divJk_hv"] = np.sum(collect_divJk_hv, axis=0) if "b" in keys_state_phys: - J_p_hv_average["J_p_v"] = np.sum(collect_J_p_v, axis=0) - J_p_hv_average["J_p_h"] = np.sum(collect_J_p_h, axis=0) + Jp_hv_average["Jp_v"] = np.sum(collect_Jp_v, axis=0) + Jp_hv_average["Jp_h"] = np.sum(collect_Jp_h, axis=0) count_final = np.sum(collect_count, axis=0) - J_k_hv_average["count"] = count_final + Jk_hv_average["count"] = count_final count_final_iso = np.sum(collect_count_iso, axis=0) - J_k_r_average["count"] = count_final_iso - S2_k_xyz_average["count2"] = count_final_iso + Jk_r_average["count"] = count_final_iso + S2_k_r_average["count2"] = count_final_iso - for index, value in np.ndenumerate(J_k_hv_average["J_k_v"]): + for index, value in np.ndenumerate(Jk_hv_average["Jk_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 + Jk_hv_average["Jk_v"][index] = 0.0 + Jk_hv_average["Jk_h"][index] = 0.0 + Jk_hv_average["Jk_l"][index] = 0.0 if "b" in keys_state_phys: - J_p_hv_average["J_p_v"][index] = 0.0 - J_p_hv_average["J_p_h"][index] = 0.0 + Jp_hv_average["Jp_v"][index] = 0.0 + Jp_hv_average["Jp_h"][index] = 0.0 else: - J_k_hv_average["J_k_v"][index] = ( + Jk_hv_average["Jk_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] + Jk_hv_average["Jk_h"][index] = ( + Jk_hv_average["Jk_h"][index] / count_final[index] + ) + Jk_hv_average["Jk_l"][index] = ( + Jk_hv_average["Jk_l"][index] / count_final[index] + ) + + Jk_hv_average["divJk_hv"][index] = ( + Jk_hv_average["divJk_hv"][index] / count_final[index] ) if "b" in keys_state_phys: - J_p_hv_average["J_p_v"][index] = ( - J_p_hv_average["J_p_v"][index] + Jp_hv_average["Jp_v"][index] = ( + Jp_hv_average["Jp_v"][index] / count_final[index] ) - J_p_hv_average["J_p_h"][index] = ( - J_p_hv_average["J_p_h"][index] + Jp_hv_average["Jp_h"][index] = ( + Jp_hv_average["Jp_h"][index] / count_final[index] ) - for index, value in np.ndenumerate(J_k_r_average["J_k_r"]): + for index, value in np.ndenumerate(Jk_r_average["Jk_average"]): if count_final_iso[index] == 0: - J_k_r_average["J_k_r"][index] = 0.0 - S2_k_xyz_average["S2_k"][index] = 0.0 + Jk_r_average["Jk_average"][index] = 0.0 + Jk_r_average["Jk_laverage"][index] = 0.0 + Jk_r_average["divJk_average"][index] = 0.0 + S2_k_r_average["S2_k_average"][index] = 0.0 else: - J_k_r_average["J_k_r"][index] = ( + Jk_r_average["Jk_average"][index] = ( value / count_final_iso[index] ) - S2_k_xyz_average["S2_k"][index] = ( - S2_k_xyz_average["S2_k"][index] + Jk_r_average["Jk_laverage"][index] = ( + Jk_r_average["Jk_laverage"][index] + / count_final_iso[index] + ) + Jk_r_average["divJk_average"][index] = ( + Jk_r_average["divJk_average"][index] + / count_final_iso[index] + ) + S2_k_r_average["S2_k_average"][index] = ( + S2_k_r_average["S2_k_average"][index] / count_final_iso[index] ) if "b" in keys_state_phys: - result = ( - J_k_r_average, - J_k_hv_average, - J_p_hv_average, - S2_k_xyz_average, - ) + result = Jk_r_average, Jk_hv_average, Jp_hv_average, S2_k_r_average + else: - result = J_k_r_average, J_k_hv_average, S2_k_xyz_average + result = Jk_r_average, Jk_hv_average, S2_k_r_average return result diff --git a/fluidsim/solvers/ns3d/output/spectra.py b/fluidsim/solvers/ns3d/output/spectra.py index 7970ccc0..6b864a9a 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_k1 is not None: - to_plot = coef_plot_k1 * ks_no0 ** (-2) * coef_norm + if coef_plot_k2 is not None: + to_plot = coef_plot_k2 * ks_no0 ** (-2) * coef_norm ax.plot(ks[1:], to_plot[1:], "k--") if xlim is not None: