diff --git a/pyTurbulence/solver.py b/pyTurbulence/solver.py
index bc1f9f70f1160e0e5476420f39b97ed2bc982ccc..d8cdbb414f9357f27aae029d562912f044d8585c 100644
--- a/pyTurbulence/solver.py
+++ b/pyTurbulence/solver.py
@@ -101,20 +101,25 @@ def solve(user_params)->dict:
domain_size[0],
domain_size[1], domain_size[2])
- _, _, tke_spectrum_dilatational = compute_tke_spectrum(ud, vd, wd,
+ _, _, tke_spectrum_dilatational = compute_tke_spectrum(u-us, v-vs, w-ws,
+ domain_size[0],
+ domain_size[1], domain_size[2])
+ _, _, tke_spectrum_total = compute_tke_spectrum(u, v, w,
domain_size[0],
domain_size[1], domain_size[2])
real_spectrum = energy_spectrum(spectrum, wave_numbers, urms, k0)
- solenoidal_spectrum_path = os.path.join(results_dir, "spectrum_solenoidal.pdf")
- dilatational_spectrum_path = os.path.join(results_dir, "spectrum_dilatational.pdf")
- plot_spectrum(wmax, wave_numbers, tke_spectrum_solenoidal, real_spectrum=real_spectrum, save_path=solenoidal_spectrum_path)
- plot_spectrum(wmax, wave_numbers, tke_spectrum_dilatational, real_spectrum=real_spectrum, save_path=dilatational_spectrum_path)
+ spectrum_path = os.path.join(results_dir, "spectrum.pdf")
+ tke_spectrum = [tke_spectrum_solenoidal, tke_spectrum_dilatational, tke_spectrum_total, real_spectrum]
+ plot_spectrum(wmax, wave_numbers, tke_spectrum, save_path=spectrum_path)
end = time.time()
spectrum_time = end - start
# Save the spectrum data
spectrum_data_path = os.path.join(results_dir, "spectrum.dat")
- np.savetxt(spectrum_data_path, np.column_stack((wave_numbers, tke_spectrum_solenoidal)), header="wave_numbers tke_spectrum")
+ np.savetxt(spectrum_data_path, np.column_stack((wave_numbers, tke_spectrum_solenoidal,
+ tke_spectrum_dilatational,
+ tke_spectrum_total)),
+ header="wave_numbers solenoidal dilatational total")
# Save the data
start = time.time()
@@ -148,8 +153,7 @@ def solve(user_params)->dict:
"pressure": pressure,
"temperature": temperature,
"wave_numbers": wave_numbers,
- "tke_spectrum_solenoidal": tke_spectrum_solenoidal,
- "tke_spectrum_dilatational": tke_spectrum_dilatational,
+ "tke_spectrum": tke_spectrum,
"real_spectrum": real_spectrum
}
diff --git a/pyTurbulence/spectrum.py b/pyTurbulence/spectrum.py
index d2c1e80482a0638796d453b7d78f450e44129e2c..a1f2406eba1a674bc52aca8c6e6e621639fe805f 100644
--- a/pyTurbulence/spectrum.py
+++ b/pyTurbulence/spectrum.py
@@ -96,16 +96,15 @@ def compute_tke_spectrum(u : np.ndarray,v : np.ndarray,w : np.ndarray,
return knyquist, kvals, tke_spectrum
-def plot_spectrum(wmax, wave_numbers, tke_spectrum, real_spectrum=None, save_path=None):
+def plot_spectrum(wmax, wave_numbers, tke_spectrum, save_path=None):
plt.figure()
- plt.loglog(wave_numbers, tke_spectrum, 'o-', label='Computed Spectrum')
- plt.axvline(x=wmax, color='r', linestyle='--', label='Cut-off Wavenumber')
-
- if real_spectrum is not None:
- plt.loglog(wave_numbers, real_spectrum, 'k', label='Real Spectrum')
-
- plt.xlabel('Wavenumber (k)')
- plt.ylabel('Energy Spectrum (E(k))')
+ plt.loglog(wave_numbers, tke_spectrum[0], 'o-', label='Solenoidal spectrum')
+ plt.loglog(wave_numbers, tke_spectrum[1], 'o-', label='Dilatational spectrum')
+ plt.loglog(wave_numbers, tke_spectrum[2], 'o-', label='Total spectrum')
+ plt.loglog(wave_numbers, tke_spectrum[3], 'k', label='Real Spectrum')
+ plt.axvline(x=wmax, color='r', linestyle='--', label='Cut-off Wavenumber')
+ plt.xlabel('k')
+ plt.ylabel('E(k)')
# plt.ylim(np.min(tke_spectrum), 1.5*np.max(tke_spectrum))
plt.ylim(1e-6, 1.5*np.max(tke_spectrum))
plt.legend()
diff --git a/pyTurbulence/syntheticTurbulence.py b/pyTurbulence/syntheticTurbulence.py
index a6afbdf35392ea5e0dd1b5d0b7551b0009dfbdf6..6a466ebab920534f3c9c952286b8c80146e5614a 100644
--- a/pyTurbulence/syntheticTurbulence.py
+++ b/pyTurbulence/syntheticTurbulence.py
@@ -168,7 +168,7 @@ def compute_solenoidal_pressure(u: np.ndarray,v: np.ndarray,w: np.ndarray,
return pressure
def compute_thermodynamic_fields(mean_density : float, mean_pressure : float, mean_temperature : float,
- u: np.ndarray, v: np.ndarray, w: np.ndarray,
+ us: np.ndarray, vs: np.ndarray, ws: np.ndarray,
ud: np.ndarray, vd: np.ndarray, wd: np.ndarray,
incompressible_pressure_fluctuations: np.ndarray,
gamma: float, Mt : float, case: int)->Tuple[np.ndarray,np.ndarray,np.ndarray,
@@ -186,11 +186,11 @@ def compute_thermodynamic_fields(mean_density : float, mean_pressure : float, me
Mean pressure.
mean_temperature : float
Mean temperature.
- u : np.ndarray
+ us : np.ndarray
Solenoidal velocity component in the x-direction.
- v : np.ndarray
+ vs : np.ndarray
Solenoidal velocity component in the y-direction.
- w : np.ndarray
+ ws : np.ndarray
Solenoidal velocity component in the z-direction.
ud : np.ndarray
Dilatational velocity component in the x-direction.
@@ -227,9 +227,9 @@ def compute_thermodynamic_fields(mean_density : float, mean_pressure : float, me
if case == 0: Mt =0.0 # no fluctuations case
compressible_pressure_fluctuations = incompressible_pressure_fluctuations * gamma * Mt**2
- u += ud*gamma*Mt**2
- v += vd*gamma*Mt**2
- w += wd*gamma*Mt**2
+ u = us + ud*gamma*Mt**2
+ v = vs + vd*gamma*Mt**2
+ w = ws + wd*gamma*Mt**2
# Assuming R = 1 and rho' = overGamma*p'
compressible_density_fluctuations = overGamma*compressible_pressure_fluctuations