Split interface and solver script

.gitignore

@@ -3,6 +3,7 @@ __pycache__/
 *.pkl
 dist/
 build/
+results/
 *.dat
 *.egg-info/
 *.egg

example.py

 import numpy as np
-from pyTurbulence.syntheticTurbulence import *
-from pyTurbulence.plot3D import *
-from pyTurbulence.spectrum import *
-import time
+from pyTurbulence.solver import solve
+from datetime import datetime
+import os
 
-params = {
-    "nbModes"           : 512, 
-    "Mt"                : 0.23,
-    "Pressure"          : 1.0,   
-    "Temperature"       : 1.0,
-    "k0"                : 8.0,
-    "domain_size"       : (2. * np.pi, 2. * np.pi, 2. * np.pi),
-    "domain_resolution" : (256,256,256),
-    "seed"              : 42,
-    "gamma"             : 1.4,
-    "case"              : 1,
-    "Spectrum"          : "PassotPouquet"
-}
+def main():
+    # Metadata
+    code_name = "~~~ pyTurbulence: Synthetic Turbulence Generator ~~~"
+    author = "Amaury Bilocq"
+    date = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
 
-# read the params dictionnary 
-nbModes = params["nbModes"]
-Mt = params["Mt"]
-mean_pressure = params["Pressure"]
-mean_temperature = params["Temperature"]
-mean_density = mean_pressure / mean_temperature
-gamma = params["gamma"]
-celerity = np.sqrt(gamma * mean_pressure / mean_density)
-Ek = 0.5*(Mt*celerity)**2
-urms = np.sqrt(2. * Ek / 3.)
-k0 = params["k0"]
-domain_size = params["domain_size"]
-domain_resolution = params["domain_resolution"]
-xi = gamma*Mt**2
-# xi = 0.0
-spectrum = params["Spectrum"]
+    # Parameters
+    params = {
+        "nbModes": 256,
+        "Mt": 0.23,
+        "Pressure": 1.0,
+        "Temperature": 1.0,
+        "k0": 8.0,
+        "domain_size": (2. * np.pi, 2. * np.pi, 2. * np.pi),
+        "domain_resolution": (64, 64, 64),
+        "seed": 42,
+        "gamma": 1.4,
+        "case": 1,
+        "Spectrum": "PassotPouquet"
+    }
 
-# Generate the solenoidal velocity field
-start = time.time()
-u,v,w,wmax = compute_solenoidal_velocities(spectrum, nbModes, urms, k0, domain_size, domain_resolution, seed=42)
-TKE = np.mean(0.5*(u.reshape(-1)**2 + v.reshape(-1)**2 + w.reshape(-1)**2))
-end = time.time()
-print(f"done computing the solenoidal velocity field with TKE = {TKE}, expected TKE = {Ek} in {end-start} seconds")
+    # Create results directory
+    results_dir = "results"
+    if not os.path.exists(results_dir):
+        os.makedirs(results_dir)
 
-# Generate the incompressible pressure fluctuations
-start = time.time()
-incompressible_pressure_fluctuations = compute_solenoidal_pressure(u,v,w, domain_resolution, domain_size)
-end = time.time()
-print(f"done computing the incompressible pressure fluctuations in {end-start} seconds")
+    # Run the solver
+    results = solve(params, results_dir)
 
-# 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"])
-end = time.time()
-print(f"done computing the thermodynamic fields in {end-start} seconds")
+    # Log file
+    log_file = os.path.join(results_dir, "log.txt")
+    with open(log_file, "w") as f:
+        # Display metadata
+        f.write(f"{'='*70}\n")
+        f.write(f"{code_name:^70}\n")
+        f.write(f"{'='*70}\n")
+        f.write(f"{'Author:':<20} {author}\n")
+        f.write(f"{'Date:':<20} {date}\n")
+        f.write(f"{'='*70}\n")
 
-# Plot 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])
-plot_spectrum(wmax, wave_numbers, tke_spectrum)
-end = time.time()
-print(f"done plotting the spectrum in {end-start} seconds")
+        # Display parameters
+        f.write(f"{'Parameters':<20}\n")
+        f.write(f"{'-'*70}\n")
+        for key, value in params.items():
+            f.write(f"{key:<20} : {value}\n")
+        f.write(f"{'='*70}\n")
 
-# Generate the visualisation 
-# plot3D(u, domain_size, domain_resolution, "U")
-# plot3D(pressure, domain_size, domain_resolution, "Pressure")
+        # Display results
+        f.write(f"{'Task':<50} {'Time (seconds)':>20}\n")
+        f.write(f"{'-'*70}\n")
+        f.write(f"{'Computing the solenoidal velocity field':<50} {results['solenoidal_time']:>20.4f}\n")
+        f.write(f"{'Computing the incompressible pressure fluctuations':<50} {results['pressure_time']:>20.4f}\n")
+        f.write(f"{'Computing the thermodynamic fields':<50} {results['thermodynamic_time']:>20.4f}\n")
+        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"{'='*70}\n")
 
-# Save the data
-start = time.time()
-u.tofile("u.dat")
-v.tofile("v.dat")
-w.tofile("w.dat")
-pressure.tofile("p.dat")
-temperature.tofile("T.dat")
-end = time.time()
-print(f"done saving the data in {end-start} seconds")
+    # Print log to console
+    with open(log_file, "r") as f:
+        print(f.read())
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file

pyTurbulence/solver.py

+import numpy as np
+from pyTurbulence.syntheticTurbulence import compute_solenoidal_velocities, compute_solenoidal_pressure, compute_thermodynamic_fields
+from pyTurbulence.spectrum import compute_tke_spectrum, plot_spectrum, energy_spectrum
+import time
+import os
+
+def solve(params, results_dir):
+    # Read the params dictionary
+    nbModes = params["nbModes"]
+    Mt = params["Mt"]
+    mean_pressure = params["Pressure"]
+    mean_temperature = params["Temperature"]
+    mean_density = mean_pressure / mean_temperature
+    gamma = params["gamma"]
+    celerity = np.sqrt(gamma * mean_pressure / mean_density)
+    Ek = 0.5 * (Mt * celerity) ** 2
+    urms = np.sqrt(2. * Ek / 3.)
+    k0 = params["k0"]
+    domain_size = params["domain_size"]
+    domain_resolution = params["domain_resolution"]
+    xi = gamma * Mt ** 2
+    spectrum = params["Spectrum"]
+
+    # Generate the solenoidal velocity field
+    start = time.time()
+    u, v, w, wmax = compute_solenoidal_velocities(spectrum, nbModes, urms, k0, domain_size, domain_resolution, seed=params["seed"])
+    TKE = np.mean(0.5 * (u.reshape(-1) ** 2 + v.reshape(-1) ** 2 + w.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, domain_resolution, domain_size)
+    end = time.time()
+    pressure_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"])
+    end = time.time()
+    thermodynamic_time = end - start
+
+    # Plot 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])
+    real_spectrum = energy_spectrum(spectrum, wave_numbers, urms, k0)
+    spectrum_path = os.path.join(results_dir, "spectrum.png")
+    plot_spectrum(wmax, wave_numbers, tke_spectrum, real_spectrum=real_spectrum, save_path=spectrum_path)
+    end = time.time()
+    spectrum_time = end - start
+
+    # Save the data
+    start = time.time()
+    u.tofile(os.path.join(results_dir, "u.dat"))
+    v.tofile(os.path.join(results_dir, "v.dat"))
+    w.tofile(os.path.join(results_dir, "w.dat"))
+    pressure.tofile(os.path.join(results_dir, "p.dat"))
+    temperature.tofile(os.path.join(results_dir, "T.dat"))
+    end = time.time()
+    save_time = end - start
+
+    return {
+        "TKE": TKE,
+        "Ek": Ek,
+        "solenoidal_time": solenoidal_time,
+        "pressure_time": pressure_time,
+        "thermodynamic_time": thermodynamic_time,
+        "spectrum_time": spectrum_time,
+        "save_time": save_time
+    }
\ No newline at end of file

pyTurbulence/spectrum.py

@@ -92,11 +92,20 @@ def compute_tke_spectrum(u,v,w,lx,ly,lz):
 	
     return knyquist, kvals, Abins/(kbins[1:] - kbins[:-1])
 
-def plot_spectrum(knyquist, wave_numbers, tke_spectrum):
+def plot_spectrum(wmax, wave_numbers, tke_spectrum, real_spectrum=None, save_path=None):
     plt.figure()
-    plt.loglog(wave_numbers, tke_spectrum,'o-')
-    plt.axvline(x=knyquist, color='r', linestyle='--')
-    plt.xlabel('k')
-    plt.ylabel('E(k)')
+    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.legend()
     plt.grid()
-    plt.show()
+    
+    if save_path:
+        plt.savefig(save_path)
+    else:
+        plt.show()

pyproject.toml

@@ -18,7 +18,8 @@ dependencies = [
     "numba>=0.56",      # For JIT optimization
     "pytest>=7.0",      # For testing
     "matplotlib>=3.4",  # For plotting
-    "scipy>=1.7"        # For scientific computations
+    "scipy>=1.7",       # For scientific computations
+    "typing-extensions>=3.10"  # For type hints and compatibility
 ]
 
 [project.optional-dependencies]