modify computation of rhs for pressure

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
     }
 

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