diff --git a/.gitignore b/.gitignore
index 0788b73ea4705086bc610d08edbc1604e2fcc062..8e4864f8b5d8cf5b96c5953f5aa2767c234e88b0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,6 +3,7 @@ __pycache__/
*.pkl
dist/
build/
+results/
*.dat
*.egg-info/
*.egg
diff --git a/example.py b/example.py
index f7bd6516b48f46d2843a042152b395659513c40b..703ef9f997d4a0e116642d3e226df77e4da94d88 100644
--- a/example.py
+++ b/example.py
@@ -1,83 +1,71 @@
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
diff --git a/pyTurbulence/solver.py b/pyTurbulence/solver.py
new file mode 100644
index 0000000000000000000000000000000000000000..53beea3c1d3f8889a98d808185f27476dbb57b17
--- /dev/null
+++ b/pyTurbulence/solver.py
@@ -0,0 +1,70 @@
+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
diff --git a/pyTurbulence/spectrum.py b/pyTurbulence/spectrum.py
index 79b64119267385caaf4bd383cc917795e17de11c..a9af613be1f0b24501d9758fd38921a4c15246cc 100644
--- a/pyTurbulence/spectrum.py
+++ b/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()
diff --git a/pyproject.toml b/pyproject.toml
index 6f15e7a0ccff278f93edc0d40fa6f362f65182d2..6ce76faca37e47b95f296fe21a92cabd59d45cb4 100644
--- a/pyproject.toml
+++ b/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]