diff --git a/README.md b/README.md
index 37036258543400cc15d20a4985c1de7ed96e62cb..2e6ef7c04cd109b8d945f7855be2cb3b0fad06c4 100644
--- a/README.md
+++ b/README.md
@@ -19,7 +19,8 @@ TODO:
To install the package, navigate to the root directory of the project and run:
```sh
-pip install .
+pip3 install -- upgrade pip
+pip3 install -e .
```
## Usage
diff --git a/example.py b/example.py
index 42601fae1ac03ccf883287b13db93d07c11eda39..cd3851b827977a25bae7ae49ecae6cdb132789ab 100644
--- a/example.py
+++ b/example.py
@@ -27,6 +27,9 @@ def main():
# Plot the 3D fields
plot3D(data["u"], params["domain_size"],
params["domain_resolution"], name="u")
+
+ plot3D(data["pressure"], params["domain_size"],
+ params["domain_resolution"], name="p")
if __name__ == "__main__":
main()
\ No newline at end of file
diff --git a/pyTurbulence/poissonSolver.py b/pyTurbulence/poissonSolver.py
index e678038ff39a74330d3c9f1f04dacf1438885dd8..f249b06d59ebc7d8c01be06538c07d4797c8b39b 100644
--- a/pyTurbulence/poissonSolver.py
+++ b/pyTurbulence/poissonSolver.py
@@ -100,4 +100,64 @@ def poisson_3d(f, nx, ny, nz, hx, hy, hz)->np.ndarray:
res[nx-1, :, :] = res[0, :, :]
res[:, ny-1, :] = res[:, 0, :]
res[:, :, nz-1] = res[:, :, 0]
- return res
\ No newline at end of file
+ return res
+
+def dilatational_velocities(theta, nx, ny, nz, hx, hy, hz)->np.ndarray:
+ """
+ Compute the dilatational velocity field using Fourier space projection.
+
+ Parameters
+ ----------
+ theta: np.ndarray
+ Divergence of the velocity field
+ nx: int
+ Number of points in the x-direction
+ ny: int
+ Number of points in the y-direction
+ nz: int
+ Number of points in the z-direction
+ hx: float
+ Spatial discretization step in the x-direction
+ hy: float
+ Spatial discretization step in the y-direction
+ hz: float
+ Spatial discretization step in the z-direction
+
+ Returns
+ -------
+ vd_x, vd_y, vd_z : np.ndarray
+ Dilatational velocity components in real space.
+ """
+
+ # Dimensions (using full grid size)
+ nnx, nny, nnz = nx, ny, nz
+ rnz = nnz//2+1
+
+ # Forward 3D FFT (using rfftn for real-valued input)
+ theta_hat = np.fft.rfftn(theta, s=(nnx, nny, nnz))
+
+ # Compute eigenvalues of gradient operator in Fourier space with zero mean
+ kx = np.fft.fftfreq(nnx, d=hx) * 2 * np.pi
+ ky = np.fft.fftfreq(nny, d=hy) * 2 * np.pi
+ kz = np.fft.fftfreq(nnz, d=hz)[:rnz] * 2 * np.pi
+
+ den = np.empty((nnx,nny,rnz),dtype=np.float64)
+ den[:, :, :] = kx[:, np.newaxis, np.newaxis]**2 + ky[np.newaxis, :, np.newaxis]**2 + kz[np.newaxis, np.newaxis, :]**2
+ den[0, 0, 0] = 1.
+
+ # Compute dilatational velocity components
+ vd_x_hat = -1j*theta_hat * kx[:, np.newaxis, np.newaxis] / den
+ vd_y_hat = -1j*theta_hat * ky[np.newaxis, :, np.newaxis] / den
+ vd_z_hat = -1j*theta_hat * kz[np.newaxis, np.newaxis, :] / den
+
+ # Ensure conjugate symmetry for real-valued input
+ vd_x_hat[0, 0, 0] = 0.
+ vd_y_hat[0, 0, 0] = 0.
+ vd_z_hat[0, 0, 0] = 0.
+
+ # Inverse 3D FFT
+ vd_x = np.fft.irfftn(vd_x_hat, s=(nnx, nny, nnz))
+ vd_y = np.fft.irfftn(vd_y_hat, s=(nnx, nny, nnz))
+ vd_z = np.fft.irfftn(vd_z_hat, s=(nnx, nny, nnz))
+
+ return vd_x, vd_y, vd_z
diff --git a/pyTurbulence/solver.py b/pyTurbulence/solver.py
index ee82f27f15b6f8c74884a18d5a2e6e801a4f80f2..9dbb6c685af65f4b76f6f170eb065cbe7e858775 100644
--- a/pyTurbulence/solver.py
+++ b/pyTurbulence/solver.py
@@ -1,5 +1,5 @@
import numpy as np
-from pyTurbulence.syntheticTurbulence import compute_solenoidal_velocities, compute_solenoidal_pressure, compute_thermodynamic_fields
+from pyTurbulence.syntheticTurbulence import *
from pyTurbulence.spectrum import compute_tke_spectrum, plot_spectrum, energy_spectrum
from pyTurbulence.plot3D import plot3D
from datetime import datetime
@@ -71,6 +71,15 @@ def solve(user_params)->dict:
end = time.time()
pressure_time = end - start
+ # Generate the dilatational velocity field
+ start = time.time()
+ ud, vd, wd = compute_dilatational_velocities(u, v, w, 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"])
@@ -105,6 +114,7 @@ def solve(user_params)->dict:
"Ek": Ek,
"solenoidal_time": solenoidal_time,
"pressure_time": pressure_time,
+ "dilatational_time": dilatational_time,
"thermodynamic_time": thermodynamic_time,
"spectrum_time": spectrum_time,
"save_time": save_time
@@ -284,6 +294,7 @@ def _write_log_file(params, results):
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 dilatational velocity field':<50} {results['dilatational_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")
diff --git a/pyTurbulence/syntheticTurbulence.py b/pyTurbulence/syntheticTurbulence.py
index 64f0db35a5dba00eec9bb005a90c6fcfb42148cd..b9d1150d5c5cb98d55f5a118a71170a302b0d490 100644
--- a/pyTurbulence/syntheticTurbulence.py
+++ b/pyTurbulence/syntheticTurbulence.py
@@ -1,5 +1,5 @@
import numpy as np
-from pyTurbulence.poissonSolver import poisson_3d
+from pyTurbulence.poissonSolver import poisson_3d, dilatational_velocities
from pyTurbulence.spectrum import energy_spectrum
from numba import jit, prange
from typing import Tuple
@@ -167,7 +167,7 @@ def compute_solenoidal_pressure(u: np.ndarray,v: np.ndarray,w: np.ndarray,
rhs = np.zeros([nx,ny,nz])
for i in range(3):
for j in range(3):
- rhs = -1.*grad[i][j]*grad[j][i]
+ rhs += -1.*grad[i][j]*grad[j][i]
pressure = poisson_3d(rhs, nx, ny, nz, hx, hy, hz)
return pressure
@@ -218,4 +218,90 @@ def compute_thermodynamic_fields(mean_density, mean_pressure, mean_temperature,
temperature = mean_temperature*(1.+compressible_temperature_fluctuations)
pressure = mean_pressure*(1.+compressible_pressure_fluctuations)
- return density, pressure, temperature
\ No newline at end of file
+ return density, pressure, temperature
+
+def compute_dilatational_velocities(u: np.ndarray,v: np.ndarray,w: np.ndarray,p: np.ndarray,
+ gamma: float, domain_resolution: Tuple[int,int,int],
+ domain_size: Tuple[float,float,float]) -> np.ndarray:
+ """
+ Compute the dilatational velocity field from the velocity components and the pressure field.
+ First, Eq.7 from Ristorcelli and Blaisdell is solved to get the dilatation field.
+ Then, the dilatational velocity field is computed using the dilatation field (Eq.12).
+
+ Parameters
+ ----------
+ 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.
+ p : np.ndarray
+ Pressure field.
+ gamma : float
+ Heat capacity ratio.
+ domain_resolution : tuple
+ Number of grid points in each direction (nx, ny, nz).
+ domain_size : tuple
+ Physical size of the domain in each direction (Lx, Ly, Lz).
+
+ Returns
+ -------
+ dilatational_velocity : np.ndarray
+ Dilatational velocity field.
+ """
+
+ nx = domain_resolution[0]
+ ny = domain_resolution[1]
+ nz = domain_resolution[2]
+
+ hx = domain_size[0]/nx
+ hy = domain_size[1]/ny
+ hz = domain_size[2]/nz
+
+ velocity = np.array([u,v,w])
+ dudx,dudy,dudz = np.gradient(u,hx,hy,hz,edge_order=2)
+ dvdx,dvdy,dvdz = np.gradient(v,hx,hy,hz,edge_order=2)
+ dwdx,dwdy,dwdz = np.gradient(w,hx,hy,hz,edge_order=2)
+ gradVelocity = np.array([[dudx,dudy,dudz],
+ [dvdx,dvdy,dvdz],
+ [dwdx,dwdy,dwdz]])
+ dpdx,dpdy,dpdz = np.gradient(p,hx,hy,hz,edge_order=2)
+ gradPressure = np.array([dpdx,dpdy,dpdz])
+
+ # Compute time derivative of pressure fluctuations (Eq.11)
+ #rhs = 2*((u_k * du_i/dx_k + dp/dx_i) * u_j)
+ rhs = np.zeros([nx,ny,nz])
+ for i in range(3):
+ for j in range(3):
+ for k in range(3):
+ rhs += 2.*(velocity[k]*gradVelocity[i][k]+gradPressure[i])*velocity[j]
+
+ # compute second derivative of rhs
+ # (p,t),jj = (rhs),ij
+ drhsdx,drhsdy,drhsdz = np.gradient(rhs,hx,hy,hz,edge_order=2)
+ drhsdxdx, drhsdxdy, drhsdxdz = np.gradient(drhsdx,hx,hy,hz,edge_order=2)
+ drhsdydx, drhsdydy, drhsdydz = np.gradient(drhsdy,hx,hy,hz,edge_order=2)
+ drhsdzdx, drhsdzdy, drhsdzdz = np.gradient(drhsdz,hx,hy,hz,edge_order=2)
+
+ hessian_rhs = np.array([[drhsdxdx,drhsdxdy,drhsdxdz],
+ [drhsdydx,drhsdydy,drhsdydz],
+ [drhsdzdx,drhsdzdy,drhsdzdz]])
+ rhs_ij = np.zeros([nx,ny,nz])
+ for i in range(3):
+ for j in range(3):
+ rhs_ij += hessian_rhs[i][j]
+
+ ddt_p = poisson_3d(rhs_ij, nx, ny, nz, hx, hy, hz)
+
+ # Compute dilatation from Eq.7
+ # -gamma*dilatation = ddt_p+v_k*dpdxk
+ dilatation = np.zeros_like(p)
+ for i in range(3):
+ dilatation += velocity[i]*gradPressure[i]
+ dilatation += ddt_p
+ dilatation /= -gamma
+
+ # Compute dilatational velocity from Eq.12
+ vd_x, vd_y, vd_z = dilatational_velocities(dilatation, nx, ny, nz, hx, hy, hz)
+ return vd_x, vd_y, vd_z
diff --git a/pyproject.toml b/pyproject.toml
index 6ce76faca37e47b95f296fe21a92cabd59d45cb4..330bd32720ae375369864094a65f060480ef61b2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,11 +7,11 @@ name = "pyTurbulence"
version = "0.1.0"
description = "Generator of solenoidal velocity and pressure fields for turbulent flows"
authors = [
- {name = "Amaury", email = "amaury.bilocq@uliege.be"}
+ {name = "Amaury Bilocq", email = "amaury.bilocq@uliege.be"}
]
license = {text = "MIT"}
readme = "README.md"
-requires-python = ">=3.6"
+requires-python = ">=3.8"
dependencies = [
"numpy>=1.21", # Required for numerical computations
@@ -24,7 +24,14 @@ dependencies = [
[project.optional-dependencies]
dev = [
- "pytest-cov" # Adds code coverage support for pytest
+ "pytest-cov", # Adds code coverage support for pytest
+ "black", # Code formatter
+ "flake8", # Linter for better code quality
+ "mypy" # Type checking
+]
+test = [
+ "pytest",
+ "pytest-cov"
]
[tool.pytest.ini_options]
@@ -32,4 +39,5 @@ addopts = "--maxfail=0 --disable-warnings"
testpaths = ["tests"]
[tool.setuptools]
-packages = ["pyTurbulence"]
\ No newline at end of file
+packages = ["pyTurbulence"]
+include-package-data = true # Ensure all package data is included
\ No newline at end of file