|
1 | 1 | import matplotlib.pyplot as plt
|
2 | 2 | import numpy as np
|
| 3 | +from colorama import Fore |
| 4 | +from pretty_print import arr_debug, print_header, print_param, print_section |
3 | 5 | from scipy.fft import fftfreq, ifft2
|
4 | 6 |
|
| 7 | +""" |
| 8 | +TRACKING: |
| 9 | +- [ ] Check if dy is correct; we are integrating in fourier space |
| 10 | +- [ ] Implement eq16 again |
| 11 | +- [ ] Set up the comparison numerical integration plot |
| 12 | +- [ ] Implement 10 realizations and plot average for F11, F22 |
| 13 | +
|
| 14 | +
|
| 15 | +For N1 = 2**10, N2 = 2**7, |
| 16 | + Field min is |
| 17 | + F11 values are |
| 18 | +""" |
| 19 | + |
| 20 | +################################ |
| 21 | +# Flags |
| 22 | +plot_field = True |
| 23 | +plot_spectra = True |
| 24 | +use_eq15 = True # If False, use eq16 |
| 25 | +large_domain = True |
| 26 | + |
5 | 27 | # Physical params
|
6 | 28 | sigma2 = 2.0
|
7 | 29 | L_2d = 15000.0
|
|
11 | 33 | c = (8 * sigma2) / (9 * L_2d ** (2 / 3))
|
12 | 34 |
|
13 | 35 | # Domain params
|
14 |
| -L1 = 40 * L_2d |
15 |
| -L2 = 5 * L_2d |
| 36 | +L1 = 40 * L_2d if large_domain else L_2d |
| 37 | +L2 = 5 * L_2d if large_domain else L_2d / 8 |
16 | 38 | N1 = 2**10
|
17 | 39 | N2 = 2**7
|
18 | 40 |
|
19 | 41 | dx = L1 / N1
|
20 | 42 | dy = L2 / N2
|
21 | 43 |
|
22 |
| - |
23 |
| -################################ |
24 |
| -# Flags |
25 |
| -plot_field = True |
26 |
| -plot_spectra = True |
27 |
| - |
28 |
| -use_eq15 = True # If False, use eq16 |
29 |
| - |
30 |
| - |
31 |
| -print("Physical parameters:") |
32 |
| -print(f"\t sigma2 = {sigma2:.2f}") |
33 |
| -print(f"\t L_2d = {L_2d:.2f}") |
34 |
| -print(f"\t psi = {psi:.2f}") |
35 |
| -print(f"\t z_i = {z_i:.2f}") |
36 |
| - |
37 |
| -print(f"Obtained c: {c:.2f}") |
38 |
| - |
39 |
| -print("Domain parameters:") |
40 |
| -print(f"\t L1 = {L1:.2f}") |
41 |
| -print(f"\t L2 = {L2:.2f}") |
42 |
| -print(f"\t N1 = {N1}") |
43 |
| -print(f"\t N2 = {N2}") |
44 |
| - |
45 |
| -print("Problem parameters:") |
46 |
| -print(f"\t Use Equation 15 = {use_eq15}") |
47 |
| -print(f"\t Plot field = {plot_field}") |
48 |
| -print(f"\t Plot spectra = {plot_spectra}") |
49 |
| - |
| 44 | +# Replace your print statements with these prettier versions |
| 45 | +print_header("WIND FIELD SIMULATOR") |
| 46 | + |
| 47 | +print_section("Physical Parameters") |
| 48 | +print_param("sigma2", f"{sigma2:.2f}", "m²/s²") |
| 49 | +print_param("L_2d", f"{L_2d:.2f}", "m") |
| 50 | +print_param("psi", f"{np.rad2deg(psi):.2f}", "degrees") |
| 51 | +print_param("z_i", f"{z_i:.2f}", "m") |
| 52 | +print_param("c", f"{c:.4f}") |
| 53 | + |
| 54 | +print_section("Domain Parameters") |
| 55 | +print_param("L1", f"{L1:.2f}", "m") |
| 56 | +print_param("L2", f"{L2:.2f}", "m") |
| 57 | +print_param("N1", N1) |
| 58 | +print_param("N2", N2) |
| 59 | +print_param("dx", f"{dx:.2f}", "m") |
| 60 | +print_param("dy", f"{dy:.2f}", "m") |
| 61 | + |
| 62 | +print_section("Problem Configuration") |
| 63 | +print_param("Use Equation 15", f"{Fore.GREEN if use_eq15 else Fore.RED}{use_eq15}") |
| 64 | +print_param("Plot field", f"{Fore.GREEN if plot_field else Fore.RED}{plot_field}") |
| 65 | +print_param("Plot spectra", f"{Fore.GREEN if plot_spectra else Fore.RED}{plot_spectra}") |
| 66 | +print_param("Large domain", f"{Fore.GREEN if large_domain else Fore.RED}{large_domain}") |
50 | 67 |
|
51 | 68 | ###############################################################################################################
|
52 | 69 | ###############################################################################################################
|
|
120 | 137 | # C_12[i, j] = ((2 * np.pi)**2 / (L1 * L2) * phi_12[i, j])
|
121 | 138 | C_12[i, j] = 1
|
122 | 139 |
|
123 |
| -exit() |
124 |
| - |
125 | 140 | # Generate Gaussian white noise
|
126 | 141 | eta_1 = np.random.normal(0, 1, size=(N1, N2)) + 1j * np.random.normal(0, 1, size=(N1, N2))
|
127 | 142 | eta_2 = np.random.normal(0, 1, size=(N1, N2)) + 1j * np.random.normal(0, 1, size=(N1, N2))
|
|
133 | 148 | u1 = np.real(ifft2((C_11 * eta_1) + (C_12 * eta_2))) # Longitudinal component
|
134 | 149 | u2 = np.real(ifft2((C_12 * eta_1) + (C_22 * eta_2))) # Transverse component
|
135 | 150 |
|
| 151 | +arr_debug(u1, "u1", plot_heatmap=False) |
| 152 | +arr_debug(u2, "u2", plot_heatmap=False) |
| 153 | + |
| 154 | + |
136 | 155 | # Verify total variance (should match sigma2)
|
137 | 156 | var_u1 = np.var(u1)
|
138 | 157 | var_u2 = np.var(u2)
|
139 |
| -print(f"Variance of u1: {var_u1:.8f}, Variance of u2: {var_u2:.8f}") |
140 |
| -print(f"Target variance: {sigma2:.4f}") |
| 158 | + |
| 159 | +# Later in your code, replace the variance print statements with: |
| 160 | +print_section("Variance Verification") |
| 161 | +print_param("Variance of u1", f"{var_u1:.8f}", "m²/s²") |
| 162 | +print_param("Variance of u2", f"{var_u2:.8f}", "m²/s²") |
| 163 | +print_param("Target variance", f"{sigma2:.4f}", "m²/s²") |
| 164 | +print_param("Ratio u1/target", f"{var_u1/sigma2:.4f}") |
| 165 | +print_param("Ratio u2/target", f"{var_u2/sigma2:.4f}") |
141 | 166 |
|
142 | 167 | # Now plot
|
143 | 168 | if plot_field:
|
144 |
| - plt.figure(figsize=(12, 5)) |
| 169 | + plt.figure(figsize=(10, 6), dpi=100) |
| 170 | + |
| 171 | + x_km = np.linspace(0, L1 / 1000, N1) |
| 172 | + y_km = np.linspace(0, L2 / 1000, N2) |
| 173 | + X_km, Y_km = np.meshgrid(x_km, y_km, indexing="ij") |
145 | 174 |
|
146 |
| - # Longitudinal (u) wind field - subplot (a) |
147 |
| - plt.subplot(211) # 1 row, 2 columns, 1st plot |
148 |
| - plt.imshow(u1.T, cmap="coolwarm") |
149 |
| - plt.colorbar(label="m/s") |
| 175 | + plt.subplot(211) |
| 176 | + im1 = plt.pcolormesh(X_km.T, Y_km.T, u1.T, cmap="RdBu_r", shading="auto") |
| 177 | + cbar1 = plt.colorbar(im1, label="[m s$^{-1}$]") |
150 | 178 | plt.xlabel("x [km]")
|
151 | 179 | plt.ylabel("y [km]")
|
152 | 180 | plt.title("(a) u")
|
153 | 181 |
|
154 |
| - # Transverse (v) wind field - subplot (b) |
155 |
| - plt.subplot(212) # 1 row, 2 columns, 2nd plot |
156 |
| - plt.imshow(u2.T, cmap="coolwarm") |
157 |
| - plt.colorbar(label="m/s") |
| 182 | + plt.subplot(212) |
| 183 | + im2 = plt.pcolormesh(X_km.T, Y_km.T, u2.T, cmap="RdBu_r", shading="auto") |
| 184 | + cbar2 = plt.colorbar(im2, label="[m s$^{-1}$]") |
158 | 185 | plt.xlabel("x [km]")
|
159 | 186 | plt.ylabel("y [km]")
|
160 | 187 | plt.title("(b) v")
|
161 | 188 |
|
162 | 189 | plt.tight_layout()
|
163 | 190 | plt.show()
|
164 | 191 |
|
| 192 | +exit() |
165 | 193 |
|
166 | 194 | #####################################################################################################
|
167 | 195 | #####################################################################################################
|
|
182 | 210 | positive_mask = k1_arr >= 0
|
183 | 211 | k1_positive = k1_arr[positive_mask] # Positive wavenumbers
|
184 | 212 |
|
185 |
| -print("k1_arr.shape: ", k1_arr.shape) |
186 |
| -print("k1.shape: ", k1.shape) |
187 |
| - |
188 |
| -print("u1_fft.shape: ", u1_fft.shape) |
189 |
| -print("psd_u1.shape: ", psd_u1.shape) |
190 |
| -print("F11.shape: ", F11.shape) |
191 |
| - |
192 |
| -print("k1_positive.shape: ", k1_positive.shape) |
193 |
| -print("positive_mask.shape: ", positive_mask.shape) |
194 |
| - |
195 | 213 | F11_positive = F11[positive_mask]
|
196 | 214 | F22_positive = F22[positive_mask]
|
197 | 215 |
|
|
0 commit comments