Change computation of density and temperature and add spectrum for dilatational component

pyTurbulence/solver.py

@@ -24,7 +24,7 @@ def solve(user_params)->dict:
         - "domain_resolution" (list, tuple, np.ndarray): Resolution of the domain, should be a 3D vector (default: [64, 64, 64]).
         - "Spectrum" (str): Type of spectrum for turbulence generation. Options: "PassotPouquet", "Gaussian" (default: "PassotPouquet").
         - "gamma" (float): Heat capacity ratio, typically between 1 and 2 (default: 1.4).
-        - "case" (int): Dimensionalisation of the density and temperature fluctuations (from Ristorcelli & Blaisdell). Options: 1, 2 (default: 1).
+        - "case" (int): Type of turbulent fields to generate. Options: 0, 1 (default: 1). 0: incompressible - 1: compressible.
         - "output_dir" (str): Directory to store output results (default: "results").
 
     Returns
@@ -60,16 +60,16 @@ def solve(user_params)->dict:
 
     # Generate the solenoidal velocity field
     start = time.time()
-    u, v, w, wmax = compute_solenoidal_velocities(spectrum, nbModes, urms, 
+    us, vs, ws, wmax = compute_solenoidal_velocities(spectrum, nbModes, urms, 
                                                   k0, domain_size, 
                                                   domain_resolution)
-    TKE = np.mean(0.5 * (u.reshape(-1) ** 2 + v.reshape(-1) ** 2 + w.reshape(-1) ** 2))
+    TKE = np.mean(0.5 * (us.reshape(-1) ** 2 + vs.reshape(-1) ** 2 + ws.reshape(-1) ** 2))
     end = time.time()
     solenoidal_time = end - start
 
     # Generate the incompressible pressure fluctuations
     start = time.time()
-    incompressible_pressure_fluctuations = compute_solenoidal_pressure(u, v, w, 
+    incompressible_pressure_fluctuations = compute_solenoidal_pressure(us, vs, ws, 
                                                                        domain_resolution, 
                                                                        domain_size)
     end = time.time()
@@ -77,37 +77,44 @@ def solve(user_params)->dict:
 
     # Generate the dilatational velocity field
     start = time.time()
-    ud, vd, wd = compute_dilatational_velocities(u, v, w, 
+    ud, vd, wd = compute_dilatational_velocities(us, vs, ws, 
                                                  incompressible_pressure_fluctuations, 
                                                  gamma, domain_resolution, domain_size)
-    u += ud*gamma*Mt**2
-    v += vd*gamma*Mt**2
-    w += wd*gamma*Mt**2
     end = time.time()
     dilatational_time = end - start
 
     # Compute the thermodynamic fields
     start = time.time()
-    density, pressure, temperature = compute_thermodynamic_fields(mean_density, 
-                                                                  mean_pressure, 
-                                                                  mean_temperature, 
-                                                                  incompressible_pressure_fluctuations, 
-                                                                  gamma, Mt, params["case"])
+    density, pressure, temperature, u, v, w = compute_thermodynamic_fields(mean_density, 
+                                                                           mean_pressure, 
+                                                                           mean_temperature, 
+                                                                           us,vs,ws,
+                                                                           ud,vd,wd,
+                                                                           incompressible_pressure_fluctuations, 
+                                                                           gamma, Mt, params["case"])
     end = time.time()
     thermodynamic_time = end - start
 
-    # Plot the spectrum
+    # Compute the spectrum
     start = time.time()
-    knyquist, wave_numbers, tke_spectrum = compute_tke_spectrum(u, v, w, domain_size[0], domain_size[1], domain_size[2])
+    knyquist, wave_numbers, tke_spectrum_solenoidal = compute_tke_spectrum(us, vs, ws, 
+                                                                            domain_size[0], 
+                                                                            domain_size[1], domain_size[2])
+
+    _, _,                 tke_spectrum_dilatational = compute_tke_spectrum(ud, vd, wd, 
+                                                                            domain_size[0], 
+                                                                            domain_size[1], domain_size[2])
     real_spectrum = energy_spectrum(spectrum, wave_numbers, urms, k0)
-    spectrum_path = os.path.join(results_dir, "spectrum.pdf")
-    plot_spectrum(wmax, wave_numbers, tke_spectrum, real_spectrum=real_spectrum, save_path=spectrum_path)
+    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)
     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)), header="wave_numbers tke_spectrum")
+    np.savetxt(spectrum_data_path, np.column_stack((wave_numbers, tke_spectrum_solenoidal)), header="wave_numbers tke_spectrum")
 
     # Save the data
     start = time.time()
@@ -141,7 +148,8 @@ def solve(user_params)->dict:
         "pressure": pressure,
         "temperature": temperature,
         "wave_numbers": wave_numbers,
-        "tke_spectrum": tke_spectrum,
+        "tke_spectrum_solenoidal": tke_spectrum_solenoidal,
+        "tke_spectrum_dilatational": tke_spectrum_dilatational,
         "real_spectrum": real_spectrum
     }
 
@@ -203,7 +211,7 @@ def _validate_params(params) -> dict:
     # Define allowed values for specific parameters
     allowed_values = {
         "Spectrum": {"PassotPouquet", "Gaussian"},  # Allowed spectrum types
-        "case": {1, 2},  # Allowed integer cases
+        "case": {0, 1},  # Allowed integer cases
     }
 
     # Define numeric ranges for specific parameters

pyTurbulence/spectrum.py

@@ -26,8 +26,8 @@ def energy_spectrum(spectrum: str, k: np.ndarray, urms: float, k0: float) -> np.
     E_k = np.empty_like(k)
     if spectrum == "PassotPouquet":
         E_k = (urms**2 * 16. * np.sqrt(2. / np.pi) * k**4 / k0**5) * np.exp(-2. * (k / k0)**2)
-    if spectrum == "Constant":
-        E_k = urms*k**4*np.exp(-2.*k**2/k0**2)
+    if spectrum == "Gaussian":
+        E_k = 0.00013*k**4*np.exp(-2.*k**2/k0**2)
     return E_k
 
 def compute_tke_spectrum(u : np.ndarray,v : np.ndarray,w : np.ndarray,
@@ -108,7 +108,8 @@ def plot_spectrum(wmax, wave_numbers, tke_spectrum, real_spectrum=None, save_pat
     
     plt.xlabel('Wavenumber (k)')
     plt.ylabel('Energy Spectrum (E(k))')
-    plt.ylim(np.min(tke_spectrum), 1.5*np.max(tke_spectrum))
+    # plt.ylim(np.min(tke_spectrum), 1.5*np.max(tke_spectrum))
+    plt.ylim(1e-6, 1.5*np.max(tke_spectrum))
     plt.legend()
     plt.grid()
     if save_path:

pyTurbulence/syntheticTurbulence.py

@@ -131,7 +131,7 @@ def compute_solenoidal_pressure(u: np.ndarray,v: np.ndarray,w: np.ndarray,
                                 domain_resolution: Tuple[int,int,int],
                                 domain_size: Tuple[float,float,float]) -> np.ndarray:
     """
-    Compute the solenoidal pressure field from the velocity components.
+    Compute the solenoidal pressure field from the solenoidal velocity components.
 
     Parameters
     ----------
@@ -168,10 +168,14 @@ 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, 
+                                 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]:
+                                 gamma: float, Mt : float, case: int)->Tuple[np.ndarray,np.ndarray,np.ndarray,
+                                                                             np.ndarray,np.ndarray,np.ndarray]:
     """
     Compute the density, pressure and temperature fields from the incompressible pressure fluctuations.
+    Combine solenoidal and dilatational velocity fields to get the total velocity field.
     Eq (4-5) + ansatz from Ristorcelli and Blaisdell : p' = gamma * Mt^2 * p1
     
     Parameters
@@ -182,6 +186,18 @@ def compute_thermodynamic_fields(mean_density : float, mean_pressure : float, me
         Mean pressure.
     mean_temperature : float
         Mean temperature.
+    u : np.ndarray
+        Solenoidal velocity component in the x-direction.
+    v : np.ndarray
+        Solenoidal velocity component in the y-direction.
+    w : np.ndarray
+        Solenoidal velocity component in the z-direction.
+    ud : np.ndarray
+        Dilatational velocity component in the x-direction.
+    vd : np.ndarray
+        Dilatational velocity component in the y-direction.
+    wd : np.ndarray
+        Dilatational velocity component in the z-direction.
     incompressible_pressure_fluctuations : np.ndarray
         Incompressible pressure fluctuations.
     gamma : float
@@ -199,23 +215,29 @@ def compute_thermodynamic_fields(mean_density : float, mean_pressure : float, me
         Pressure field.
     temperature : np.ndarray
         Temperature field.
+    u : np.ndarray
+        Velocity component in the x-direction.
+    v : np.ndarray
+        Velocity component in the y-direction.
+    w : np.ndarray
+        Velocity component in the z-direction
     """
-    compressible_pressure_fluctuations = incompressible_pressure_fluctuations * gamma * Mt**2 
+    overGamma = 1./gamma
 
-    if case == 1:
-        compressible_density_fluctuations = compressible_pressure_fluctuations
-        compressible_temperature_fluctuations = compressible_pressure_fluctuations
+    if case == 0: Mt =0.0 # no fluctuations case
+    compressible_pressure_fluctuations = incompressible_pressure_fluctuations * gamma * Mt**2 
 
-    if case == 2:
-        compressible_density_fluctuations = compressible_pressure_fluctuations*Mt**2
-        compressible_temperature_fluctuations = compressible_pressure_fluctuations*Mt**2*(gamma-1)
+    u += ud*gamma*Mt**2
+    v += vd*gamma*Mt**2
+    w += wd*gamma*Mt**2
 
-    # Get fields of density - temperature - pressure
-    density     = mean_density*(1.+compressible_density_fluctuations)
-    temperature = mean_temperature*(1.+compressible_temperature_fluctuations)
+    # Assuming R = 1 and rho' = overGamma*p'
+    compressible_density_fluctuations = overGamma*compressible_pressure_fluctuations
+    density     = mean_density*(1.+compressible_density_fluctuations) 
     pressure    = mean_pressure*(1.+compressible_pressure_fluctuations)
+    temperature = pressure/density
 
-    return density, pressure, temperature
+    return density, pressure, temperature, u, v, w
 
 def compute_dilatational_velocities(u: np.ndarray,v: np.ndarray,w: np.ndarray,p: np.ndarray,
                                     gamma: float, domain_resolution: Tuple[int,int,int],
@@ -277,7 +299,7 @@ def compute_dilatational_velocities(u: np.ndarray,v: np.ndarray,w: np.ndarray,p:
 
     # Compute dilatation from Eq.7
     # -gamma*dilatation = ddt_p+v_k*dpdxk
-    dilatation = -overGamma*(np.sum(velocity * gradPressure, axis=0) + ddt_p)
+    dilatation = -overGamma*(np.einsum('k...,k...->...',velocity,gradPressure) + ddt_p)
 
     # Compute dilatational velocity from Eq.12
     vd_x, vd_y, vd_z = dilatational_velocities(dilatation, nx, ny, nz, hx, hy, hz)