Split TKE contributions

pyTurbulence/solver.py

@@ -55,15 +55,15 @@ def solve(user_params)->dict:
     # Compute additional parameters
     mean_density = mean_pressure / mean_temperature
     celerity = np.sqrt(gamma * mean_pressure / mean_density)
-    Ek = 0.5 * (Mt * celerity) ** 2
-    urms = np.sqrt(2. * Ek / 3.)
+    TKE_expected = 0.5 * (Mt * celerity) ** 2
+    urms = np.sqrt(2. * TKE_expected / 3.)
 
     # Generate the solenoidal velocity field
     start = time.time()
     us, vs, ws, wmax = compute_solenoidal_velocities(spectrum, nbModes, urms, 
                                                   k0, domain_size, 
                                                   domain_resolution)
-    TKE = np.mean(0.5 * (us.reshape(-1) ** 2 + vs.reshape(-1) ** 2 + ws.reshape(-1) ** 2))
+    TKE_solenoidal = np.mean(0.5 * (us.reshape(-1) ** 2 + vs.reshape(-1) ** 2 + ws.reshape(-1) ** 2))
     end = time.time()
     solenoidal_time = end - start
 
@@ -92,6 +92,9 @@ def solve(user_params)->dict:
                                                                            ud,vd,wd,
                                                                            incompressible_pressure_fluctuations, 
                                                                            gamma, Mt, params["case"])
+    TKE_total = np.mean(0.5 * (u.reshape(-1) ** 2 + v.reshape(-1) ** 2 + w.reshape(-1) ** 2))
+    TKE_dilatational = TKE_total - TKE_solenoidal
+    TKE = [TKE_solenoidal, TKE_dilatational, TKE_total, TKE_expected]
     end = time.time()
     thermodynamic_time = end - start
 
@@ -133,7 +136,6 @@ def solve(user_params)->dict:
 
     results = {
         "TKE": TKE,
-        "Ek": Ek,
         "solenoidal_time": solenoidal_time,
         "pressure_time": pressure_time,
         "dilatational_time": dilatational_time,
@@ -225,7 +227,7 @@ def _validate_params(params) -> dict:
         "Pressure": (0, None),  # Pressure must be positive
         "Temperature": (0, None),  # Temperature must be positive
         "k0": (1, None),  # Characteristic wavenumber must be >= 1
-        "gamma": (0, None),  # Heat capacity ratio is typically between 1 and 2
+        "gamma": (0, None)  # Heat capacity ratio is typically between 1 and 2
     }
 
     validated_params = default_params.copy()
@@ -321,8 +323,10 @@ def _write_log_file(params, results):
         f.write(f"{'Plotting the spectrum':<50} {results['spectrum_time']:>20.4f}\n")
         f.write(f"{'Saving the data':<50} {results['save_time']:>20.4f}\n")
         f.write(f"{'-'*70}\n")
-        f.write(f"{'Total TKE':<50} {results['TKE']:>20.4f}\n")
-        f.write(f"{'Expected TKE':<50} {results['Ek']:>20.4f}\n")
+        f.write(f"{'Solenoidal TKE':<50} {results['TKE'][0]:>20.4f}\n")
+        f.write(f"{'Dilatational TKE':<50} {results['TKE'][1]:>20.4f}\n")
+        f.write(f"{'Total TKE':<50} {results['TKE'][2]:>20.4f}\n")
+        f.write(f"{'Expected TKE':<50} {results['TKE'][3]:>20.4f}\n")
         f.write(f"{'='*70}\n")
 
     # Print log to console

pyTurbulence/spectrum.py

@@ -60,8 +60,8 @@ def compute_tke_spectrum(u : np.ndarray,v : np.ndarray,w : np.ndarray,
 
     # --- Size
     n0, n1, n2 = u.shape
-    L           = [lx, ly, lz]
-    oneOverN = 1./(n0*n1*n2)
+    L          = [lx, ly, lz]
+    oneOverN   = 1./(n0*n1*n2)
 		
     # --- FFT
     fourier_u          = np.fft.fftn(u)*oneOverN