From f43ca35f092579935ff3df1df84bf9777496768b Mon Sep 17 00:00:00 2001
From: Amaury Bilocq <amaurybilocq@mbp-de-amaury-3.user.uliege.priv>
Date: Fri, 28 Feb 2025 15:22:48 +0100
Subject: [PATCH] modify computation of rhs for pressure
---
pyTurbulence/solver.py | 4 ++--
pyTurbulence/syntheticTurbulence.py | 13 +++----------
2 files changed, 5 insertions(+), 12 deletions(-)
diff --git a/pyTurbulence/solver.py b/pyTurbulence/solver.py
index d3131e7..250e401 100644
--- a/pyTurbulence/solver.py
+++ b/pyTurbulence/solver.py
@@ -22,7 +22,7 @@ def solve(user_params)->dict:
- "k0" (int): Peak wavenumber (default: 4).
- "domain_size" (list, tuple, np.ndarray): Size of the simulation domain, should be a 3D vector (default: [2Ï€, 2Ï€, 2Ï€]).
- "domain_resolution" (list, tuple, np.ndarray): Resolution of the domain, should be a 3D vector (default: [64, 64, 64]).
- - "Spectrum" (str): Type of spectrum for turbulence generation. Options: "PassotPouquet", "Constant" (default: "PassotPouquet").
+ - "Spectrum" (str): Type of spectrum for turbulence generation. Options: "PassotPouquet", "Gaussian" (default: "PassotPouquet").
- "gamma" (float): Heat capacity ratio, typically between 1 and 2 (default: 1.4).
- "case" (int): Dimensionalisation of the density and temperature fluctuations (from Ristorcelli & Blaisdell). Options: 1, 2 (default: 1).
- "output_dir" (str): Directory to store output results (default: "results").
@@ -202,7 +202,7 @@ def _validate_params(params) -> dict:
# Define allowed values for specific parameters
allowed_values = {
- "Spectrum": {"PassotPouquet", "Gaussian", "Exponential"}, # Allowed spectrum types
+ "Spectrum": {"PassotPouquet", "Gaussian"}, # Allowed spectrum types
"case": {1, 2}, # Allowed integer cases
}
diff --git a/pyTurbulence/syntheticTurbulence.py b/pyTurbulence/syntheticTurbulence.py
index ee09627..66a0898 100644
--- a/pyTurbulence/syntheticTurbulence.py
+++ b/pyTurbulence/syntheticTurbulence.py
@@ -160,16 +160,9 @@ def compute_solenoidal_pressure(u: np.ndarray,v: np.ndarray,w: np.ndarray,
hy = domain_size[1]/ny
hz = domain_size[2]/nz
- 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)
-
- grad = np.array([[dudx,dudy,dudz],[dvdx,dvdy,dvdz],[dwdx,dwdy,dwdz]])
-
- rhs = np.zeros([nx,ny,nz])
- for i in range(3):
- for j in range(3):
- rhs += -1.*grad[i][j]*grad[j][i]
+ velocity = np.array([u,v,w])
+ gradVelocity = np.array([np.gradient(vel, hx, hy, hz, edge_order=2) for vel in velocity])
+ rhs = -np.einsum('ij...,ji...->...', gradVelocity, gradVelocity)
pressure = poisson_3d(rhs, nx, ny, nz, hx, hy, hz)
return pressure
--
GitLab