Some modifications

README.md

@@ -9,7 +9,7 @@ TODO:
 
 ## Features
 
-- Generate solenoidal velocity fields using the method of Rogallo.
+- Generate solenoidal velocity fields using Kraichnan's method, improved for mass conservation by Saad.
 - Compute incompressible pressure fluctuations using the method of Ristorcelli & Blaisdell.
 - Compute thermodynamic fields (density, pressure, temperature) from the generated velocity and pressure fields.
 - Plot the energy spectrum of the generated turbulence.
@@ -41,8 +41,8 @@ This project is licensed under the MIT License. See the LICENSE file for details
 
 ## Acknowledgments
 
-- Ristorcelli, J. R., & Blaisdell, G. A. (1997). Consistent initial conditions for the DNS of compressible turbulence. Journal of Computational Physics, 140(2), 427-466.
-- Sagaut, P. (2018). Homogeneous Turbulence Dynamics. Springer.
-- Tony Saad, Derek Cline, Rob Stoll, and James C. Sutherland. “Scalable Tools for Generating Synthetic Isotropic Turbulence with Arbitrary Spectra”. http://dx.doi.org/10.2514/1.J055230. (available online: Aug 25, 2016).
-- Saad, T., & Sutherland, J. C. (2016). Comment on “Diffusion by a random velocity field” [Phys. Fluids 13, 22 (1970)]. Physics of Fluids (1994-Present), 28(11), 119101. https://doi.org/10.1063/1.4968528.
-- Austin Richards, Tony Saad, and James C. Sutherland. “A Fast Turbulence Generator using Graphics Processing Units”, 2018 Fluid Dynamics Conference, AIAA AVIATION Forum, (AIAA 2018-3559). https://doi.org/10.2514/6.2018-3559.
+[1] Kraichnan R. H., “Diffusion by a Random Velocity Field,” Physics of Fluids, Vol. 13, No. 1, 1970, Paper 22. doi: https://doi.org/10.1063/1.1692799
+
+[2] Tony Saad, Derek Cline, Rob Stoll, and James C. Sutherland. “Scalable Tools for Generating Synthetic Isotropic Turbulence with Arbitrary Spectra”. http://dx.doi.org/10.2514/1.J055230. (available online: Aug 25, 2016).
+
+[3] Ristorcelli, J. R., & Blaisdell, G. A. (1997). Consistent initial conditions for the DNS of compressible turbulence. Journal of Computational Physics, 140(2), 427-466.
\ No newline at end of file

pyTurbulence/solver.py

@@ -183,7 +183,7 @@ def _validate_params(params) -> dict:
     # Define allowed values for specific parameters
     allowed_values = {
         "Spectrum": {"PassotPouquet", "Gaussian", "Exponential"},  # Allowed spectrum types
-        "case": {0, 1, 2},  # Allowed integer cases
+        "case": {1, 2},  # Allowed integer cases
     }
 
     # Define numeric ranges for specific parameters

pyTurbulence/spectrum.py

 import numpy as np
 import matplotlib.pyplot as plt
-from numba import njit
 from scipy import stats
 
-@njit
 def energy_spectrum(spectrum: str, k: np.ndarray, urms: float, k0: float) -> np.ndarray:
     """
     Compute the energy spectrum E(k) for a given vector of wavenumber k.
@@ -32,7 +30,7 @@ def energy_spectrum(spectrum: str, k: np.ndarray, urms: float, k0: float) -> np.
         E_k = urms*k**4*np.exp(-2.*k**2/k0**2)
     return E_k
 
-def compute_tke_spectrum(u,v,w,lx,ly,lz):
+def compute_tke_spectrum(u,v,w,lx,ly,lz,smooth=True):
     """
     Compute the turbulent kinetic energy (TKE) spectrum from the velocity fields.
 
@@ -87,10 +85,17 @@ def compute_tke_spectrum(u,v,w,lx,ly,lz):
 
     # --- Spectrum 
     Abins, _, _        = stats.binned_statistic(knrm, fourier_amplitudes, statistic = "sum", bins = kbins)
+    tke_spectrum = Abins/(kbins[1:] - kbins[:-1])
 
     knyquist = knrm * min(n0, n1, n2) / 2
+    
+    if smooth:
+        window = np.ones(int(5)) / float(5)
+        tkespecsmooth =  np.convolve(tke_spectrum, window, 'same')
+        tkespecsmooth[0:4] = tke_spectrum[0:4]  # get the first 4 values from the original data
+        tke_spectrum = tkespecsmooth
 	
-    return knyquist, kvals, Abins/(kbins[1:] - kbins[:-1])
+    return knyquist, kvals, tke_spectrum
 
 def plot_spectrum(wmax, wave_numbers, tke_spectrum, real_spectrum=None, save_path=None):
     plt.figure()

pyTurbulence/syntheticTurbulence.py

@@ -21,14 +21,14 @@ def _compute_velocity_field(u_, v_, w_, nx, ny, nz, kx, ky, kz, xc, yc, zc, um,
 
 # MIT License
 # Original function from turboGenPy by Tony Saad
-# Modified by Amaury Bilocq in 2024
+# Modified by Amaury Bilocq in 2025
 def compute_solenoidal_velocities(spectrum: str, nbModes: int, urms: float,
                                     k0: float, domain_size: Tuple[float, float, float], 
                                     domain_resolution: Tuple[int, int, int]
                                     ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, float]:
 
     """
-    Generate a synthetic turbulence velocity field using the Tony Saad methodology.
+    Generate a synthetic turbulence velocity field using the Tony Saad's methodology.
 
     Parameters
     ----------
@@ -64,12 +64,12 @@ def compute_solenoidal_velocities(spectrum: str, nbModes: int, urms: float,
     dz = Lz/nz
 
     # Proper wavenumber scaling
-    wmin = 2.*np.pi/Lx
-    wmax = np.pi/dx
+    # unique_kk, wave_instances, wave_width = wavespace(nbModes)
+    wmin = 2.*np.pi/np.min([Lx,Ly,Lz])
+    wmax = np.pi/np.min([dx,dy,dz])
 
     dk = (wmax - wmin)/nbModes
     wn = wmin + np.arange(0,nbModes)*dk
-
     dkn = np.ones(nbModes)*dk
 
     # Generate random angles - phases
@@ -109,7 +109,7 @@ def compute_solenoidal_velocities(spectrum: str, nbModes: int, urms: float,
     Ek = energy_spectrum(spectrum, km, urms, k0)
     Ek = Ek.clip(0.0)
 
-    # generate turbulence at cell centers
+    # generate turbulence
     um = np.sqrt(Ek * dkn)
     u_ = np.zeros((nx, ny, nz))
     v_ = np.zeros((nx, ny, nz))
@@ -149,6 +149,7 @@ def compute_solenoidal_pressure(u: np.ndarray,v: np.ndarray,w: np.ndarray,
     pressure : np.ndarray
         Solenoidal pressure field.
     """
+
     nx = domain_resolution[0]
     ny = domain_resolution[1]
     nz = domain_resolution[2]