First test computation of dilatational velocity field

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

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

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

pyTurbulence/solver.py

 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")

pyTurbulence/syntheticTurbulence.py

 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

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