diff --git a/pyTurbulence/solver.py b/pyTurbulence/solver.py index d3131e7395f15ecac937dd1e782f3109023cbdd7..250e4012e211b88bac6ba8c0338f4c43449d2047 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 ee0962739212f2063e5ed352280a9176cb5a579a..66a089864666c5def57f7188999aa286afca2e96 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