Merge the dilatational/spectrum/total spectrum together

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
     }
 

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

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