Skip to content

Commit e467e7d

Browse files
mdmeekermdmeeker
committed
CHECK DIFFS HERE IF TESTS FAIL
1 parent a5314e0 commit e467e7d

8 files changed

+146
-139
lines changed

‎drdmannturb/common.py

+17-9
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
"""
2-
This module contains implementations of common functions, specifically the Mann eddy lifetime function and the von Karman energy spectrum.
2+
This module contains implementations of common functions, specifically the Mann eddy lifetime function
3+
and the von Karman energy spectrum.
34
"""
45

56
__all__ = ["VKEnergySpectrum", "MannEddyLifetime", "Mann_linear_exponential_approx"]
@@ -43,12 +44,14 @@ def MannEddyLifetime(kL: Union[torch.Tensor, np.ndarray]) -> torch.Tensor:
4344
Implementation of the full Mann eddy lifetime function, of the form
4445
4546
.. math::
46-
\tau^{\mathrm{IEC}}(k)=\frac{(k L)^{-\frac{2}{3}}}{\sqrt{{ }_2 F_1\left(1 / 3,17 / 6 ; 4 / 3 ;-(k L)^{-2}\right)}}
47+
\tau^{\mathrm{IEC}}(k)=\frac{(k L)^{-\frac{2}{3}}}{\sqrt{{ }_2 F_1\left(1/3, 17/6; 4/3 ;-(kL)^{-2}\right)}}
4748
4849
This function can execute with input data that are either in Torch or numpy. However,
4950
5051
.. warning::
51-
This function depends on SciPy for evaluating the hypergeometric function, meaning a GPU tensor will be returned to the CPU for a single evaluation and then converted back to a GPU tensor. This incurs a substantial loss of performance.
52+
This function depends on SciPy for evaluating the hypergeometric function, meaning a GPU tensor will be returned
53+
to the CPU for a single evaluation and then converted back to a GPU tensor. This incurs a substantial loss of
54+
performance.
5255
5356
Parameters
5457
----------
@@ -76,12 +79,17 @@ def Mann_linear_exponential_approx(
7679
.. math::
7780
\frac{x^{-\frac{2}{3}}}{\sqrt{{ }_2 F_1\left(1 / 3,17 / 6 ; 4 / 3 ;-x^{-2}\right)}}
7881
79-
via an exponential function, which is a reasonable approximation since the resulting :math:`\tau` is nearly linear on a log-log plot. The resulting approximation is the function
82+
via an exponential function, which is a reasonable approximation since the resulting :math:`\tau` is nearly linear
83+
on a log-log plot. The resulting approximation is the function
8084
8185
.. math::
8286
\exp \left( \alpha kL + \beta \right)
8387
84-
where :math:`\alpha, \beta` are obtained from a linear regression on the hypergeometric function on the domain of interest. In particular, using this function requires that a linear regression has already been performed on the basis of the above function depending on the hypergeometric function, which is an operation performed once on the CPU. The rest of this subroutine is on the GPU and unlike the full hypergeometric approximation, will not incur any slow-down of the rest of the spectra fitting.
88+
where :math:`\alpha, \beta` are obtained from a linear regression on the hypergeometric function on the domain of
89+
interest. In particular, using this function requires that a linear regression has already been performed on the
90+
basis of the above function depending on the hypergeometric function, which is an operation performed once on the
91+
CPU. The rest of this subroutine is on the GPU and unlike the full hypergeometric approximation, will not incur
92+
any slow-down of the rest of the spectra fitting.
8593
8694
Parameters
8795
----------
@@ -112,7 +120,9 @@ def find_class(self, module, name):
112120

113121

114122
def plot_loss_logs(log_file: Union[str, Path]):
115-
"""Returns a full plot of all loss terms for a specific training log generated by TensorBoard. This is an auxiliary method and should only be used for quick visualization of the training process, the suggested method for visualizing this information is through the TensorBoard API.
123+
"""Returns a full plot of all loss terms for a specific training log generated by TensorBoard. This is an auxiliary
124+
method and should only be used for quick visualization of the training process, the suggested method for visualizing
125+
this information is through the TensorBoard API.
116126
117127
Parameters
118128
----------
@@ -127,9 +137,7 @@ def plot_loss_logs(log_file: Union[str, Path]):
127137

128138
vals_tot = {}
129139
for t_scalar in training_scalars:
130-
_, _, vals_curr = zip(
131-
*[astuple(event) for event in event_acc.Scalars(t_scalar)]
132-
)
140+
_, _, vals_curr = zip(*[astuple(event) for event in event_acc.Scalars(t_scalar)])
133141

134142
vals_tot[t_scalar] = vals_curr[1:]
135143

‎drdmannturb/fluctuation_generation/fluctuation_field_generator.py

+44-49
Original file line numberDiff line numberDiff line change
@@ -80,7 +80,6 @@ def __init__(
8080
time_scale: Optional[float] = None,
8181
energy_spectrum_scale: Optional[float] = None,
8282
path_to_parameters: Optional[Union[str, PathLike]] = None,
83-
lowfreq_params: Optional[dict] = None,
8483
seed: Optional[int] = None,
8584
blend_num=10,
8685
):
@@ -162,73 +161,69 @@ def __init__(
162161
L = length_scale
163162
Gamma = time_scale
164163

165-
# define margins and buffer
166-
time_buffer = 3 * Gamma * L
167-
spatial_margin = 1 * L # TODO: where is this used?
168-
169-
# TODO: this should not be needed
170-
# grid_levels = [level.GetInt() for level in grid_levels]
164+
self.grid_dimensions = grid_dimensions
165+
self.grid_levels = grid_levels
166+
self.blend_num = blend_num
167+
print("Grid levels is ", self.grid_levels)
168+
print("Grid dimensions is ", self.grid_dimensions)
169+
print("Blend number is ", self.blend_num)
171170

172171
# Expand on grid_levels parameter to get grid node counts in each direction
173172
grid_node_counts = 2**grid_levels + 1
174-
Nx, Ny, Nz = grid_node_counts
173+
print("Grid node counts is ", grid_node_counts)
175174

176-
# Use grid_dimensions and grid_node_counts to get grid spacings in each direction
177-
grid_spacings = [dim / size for dim, size in zip(grid_dimensions, grid_node_counts)]
178-
hx, hy, hz = grid_spacings
175+
# Obtain spacing between grid points, split node counts into Nx, Ny, Nz
176+
dx, dy, dz = (L_i / N_i for L_i, N_i in zip(self.grid_dimensions, grid_node_counts))
177+
Nx, Ny, Nz = grid_node_counts
178+
del grid_node_counts
179179

180180
# Calculate buffer and margin sizes
181-
n_buffer = ceil(time_buffer / hx)
182-
n_margins = [
183-
ceil(spatial_margin / hy), # y margin
184-
ceil(spatial_margin / hz), # z margin
185-
]
181+
# NOTE: buffer scale 3 * Gamma * L is arbitrary. Could/should be tunable param?
182+
self.n_buffer = ceil((3 * Gamma * L) / dx)
183+
184+
self.n_margin_y, self.n_margin_z = ceil(L / dy), ceil(L / dz)
185+
186+
## Spatial margin is just the length scale L
186187

187188
# Calculate shapes
188-
wind_shape = [0, Ny, Nz, 3]
189+
buffer_extension = 2 * self.n_buffer + (self.blend_num - 1 if self.blend_num > 0 else 0)
190+
margin_extension = [
191+
2 * self.n_margin_y,
192+
2 * self.n_margin_z,
193+
]
189194

190-
buffer_extension = 2 * n_buffer + (blend_num - 1 if blend_num > 0 else 0)
191-
margin_extension = [2 * margin for margin in n_margins]
195+
print("Buffer extension is ", buffer_extension)
196+
print("Margin extension is ", margin_extension)
192197

193-
noise_shape = [
198+
self.noise_shape = [
194199
Nx + buffer_extension,
195200
Ny + margin_extension[0],
196201
Nz + margin_extension[1],
197202
3,
198203
]
204+
self.new_part_shape = [Nx, Ny + margin_extension[0], Nz + margin_extension[1], 3]
205+
self.central_part = [
206+
slice(self.n_buffer, -self.n_buffer),
207+
slice(self.n_margin_y, -self.n_margin_y),
208+
slice(0, -2 * self.n_margin_z),
209+
slice(None),
210+
]
211+
self.new_part = [
212+
slice(2 * self.n_buffer + max(0, self.blend_num - 1), None),
213+
slice(None),
214+
slice(None),
215+
slice(None),
216+
]
199217

200-
new_part_shape = [Nx] + [size + margin_ext for size, margin_ext in zip([Ny, Nz], margin_extension)] + [3]
201-
202-
# Calculate slices for different parts
203-
def make_slice(start=None, end=None):
204-
return slice(start, end)
205-
206-
central_part = [make_slice()] * 4 # Initially all are full slices
207-
central_part[0] = make_slice(n_buffer, -n_buffer)
208-
central_part[1] = make_slice(n_margins[0], -n_margins[0])
209-
central_part[2] = make_slice(0, -2 * n_margins[1])
210-
211-
new_part = [make_slice()] * 4
212-
if blend_num > 0:
213-
new_part[0] = make_slice(2 * n_buffer + (blend_num - 1), None)
214-
else:
215-
new_part[0] = make_slice(2 * n_buffer, None)
216-
217-
self.grid_dimensions = grid_dimensions
218-
self.grid_levels = grid_levels
219-
220-
self.new_part = new_part
221218
self.Nx = Nx
222-
self.blend_num = blend_num
223-
self.central_part = central_part
224-
self.new_part_shape = new_part_shape
225-
self.noise_shape = noise_shape
226-
self.n_buffer = n_buffer
227-
self.n_marginy = n_margins[0]
228-
self.n_marginz = n_margins[1]
229219
self.seed = seed
220+
# NOTE: self.noise is a placeholder for now
221+
# - Also, total_fluctuation is the "accumulator" for generated blocks
230222
self.noise = None
231-
self.total_fluctuation = np.zeros(wind_shape)
223+
self.total_fluctuation = np.zeros([0, Ny, Nz, 3])
224+
225+
print("DONE WITH FLUCTUATION FIELD GENERATOR INIT")
226+
print("\n\n\n")
232227

233228
CovarianceType = Union[
234229
type[VonKarmanCovariance],

‎drdmannturb/fluctuation_generation/gaussian_random_fields.py

+3-5
Original file line numberDiff line numberDiff line change
@@ -92,7 +92,8 @@ def __init__(
9292
# TODO: This is a strange block of code
9393
if np.isscalar(grid_level):
9494
if not np.isscalar(grid_shape):
95-
raise ValueError("grid_level and grid_shape must have the same dimensions.")
95+
raise ValueError("grid_level is scalar, so grid_shape should be as well.")
96+
9697
h = 1 / 2**grid_level
9798
self.grid_shape = np.array([grid_shape] * ndim)
9899
else:
@@ -101,10 +102,10 @@ def __init__(
101102

102103
h = grid_dimensions / (2**grid_level + 1)
103104
self.grid_shape = np.array(grid_shape[:ndim])
105+
104106
self.L = h * self.grid_shape
105107

106108
### Extended window (NOTE: extension is done outside)
107-
# TODO: This is also a strange block of code. Vestigial?
108109
N_margin = 0
109110
self.ext_grid_shape = self.grid_shape
110111
self.nvoxels = self.ext_grid_shape.prod()
@@ -113,9 +114,6 @@ def __init__(
113114
### Covariance kernel
114115
self.Covariance = Covariance
115116

116-
### Sampling method
117-
self.setSamplingMethod(sampling_method)
118-
119117
# Pseudo-random number generator
120118
self.prng = np.random.RandomState()
121119
self.noise_std = np.sqrt(np.prod(h))

‎drdmannturb/fluctuation_generation/sampling_methods.py

+6
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,12 @@ def __init__(self, RandomField):
4747
super().__init__(RandomField)
4848
L, Nd, d = self.L, self.Nd, self.ndim
4949
self.Frequencies = [(2 * np.pi / L[j]) * (Nd[j] * fft.fftfreq(Nd[j])) for j in range(d)]
50+
51+
print("Sampling_method_freq init: L", L)
52+
print("Sampling_method_freq init: Nd", Nd)
53+
print("Sampling_method_freq init: d", d)
54+
print("Sampling_method_freq init: self.Frequencies", self.Frequencies)
55+
5056
self.TransformNorm = np.sqrt(L.prod())
5157
self.Spectrum = RandomField.Covariance.precompute_Spectrum(self.Frequencies)
5258

‎drdmannturb/fluctuation_generation/wind_plot.py

+27-20
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
"""
2-
Common utilities and Plotly integration for visualizing generated wind field.
2+
Common utilities and Plotly integration for visualizing generated wind field.
33
"""
44

55
import numpy as np
@@ -9,22 +9,24 @@
99
from .cmap_util import FIELD_COLORSCALE
1010

1111

12-
def create_grid(
13-
spacing: tuple[float, float, float], shape: tuple[int, int, int]
14-
) -> np.ndarray:
15-
"""Creates a 3D grid (meshgrid) from given spacing between grid points and desired shape (which should match the shape of the generated wind field, for example).
12+
def create_grid(spacing: tuple[float, float, float], shape: tuple[int, int, int]) -> np.ndarray:
13+
"""Creates a 3D grid (meshgrid) from given spacing between grid points and desired shape (which should match the
14+
shape of the generated wind field, for example).
1615
1716
Parameters
1817
----------
1918
spacing : tuple[float, float, float]
20-
Spacing array that determines the spacing of points to be used in each dimension of the 3D field. Typically, of the form grid_dimensions (a 3x1 vector representing the dimensions of the domain) divided by the grid_levels, which determine the resolution of the wind field in each respective dimension.
19+
Spacing array that determines the spacing of points to be used in each dimension of the 3D field. Typically, of
20+
the form grid_dimensions (a 3x1 vector representing the dimensions of the domain) divided by the grid_levels,
21+
which determine the resolution of the wind field in each respective dimension.
2122
shape : tuple[int, int, int]
2223
Number of points in each dimension.
2324
2425
Returns
2526
-------
2627
np.ndarray
27-
np.meshgrid object consisting of points at the provided spacing and with the specified counts in each dimension. This is 'ij' indexed (not Cartesian!).
28+
np.meshgrid object consisting of points at the provided spacing and with the specified counts in each dimension.
29+
This is 'ij' indexed (not Cartesian!).
2830
"""
2931
x = np.array([spacing[0] * n for n in range(shape[0])])
3032
y = np.array([spacing[1] * n for n in range(shape[1])])
@@ -58,18 +60,22 @@ def plot_velocity_components(
5860
surface_count=25,
5961
reshape=True,
6062
) -> go.Figure:
61-
"""Plots x, y, z components of given wind field over provided spacing. Note that the same spacing is used for all 3 velocity components.
63+
"""Plots x, y, z components of given wind field over provided spacing. Note that the same spacing is used for all 3
64+
velocity components.
6265
6366
Parameters
6467
----------
6568
spacing : tuple[float, float, float]
66-
Spacing array that determines the spacing of points to be used in each dimension of the 3D field. Typically, of the form grid_dimensions (a 3x1 vector representing the dimensions of the domain) divided by the grid_levels, which determine the resolution of the wind field in each respective dimension.
69+
Spacing array that determines the spacing of points to be used in each dimension of the 3D field. Typically, of
70+
the form grid_dimensions (a 3x1 vector representing the dimensions of the domain) divided by the grid_levels,
71+
which determine the resolution of the wind field in each respective dimension.
6772
wind_field : np.ndarray
6873
3D wind field, typically of shape (Nx, Ny, Nz, 3) (not C-layout, to be reshaped).
6974
surface_count : int, optional
7075
Number of surfaces to be used for each velocity component, by default 25
7176
reshape : bool, optional
72-
Whether to re-format the given wind field into C-order, typically the desirable choice to match the order of entries of the wind field and the provided spacing, by default True
77+
Whether to re-format the given wind field into C-order, typically the desirable choice to match the order of
78+
entries of the wind field and the provided spacing, by default True
7379
7480
Returns
7581
-------
@@ -167,20 +173,25 @@ def plot_velocity_magnitude(
167173
reshape=True,
168174
transparent=False,
169175
) -> go.Figure:
170-
"""Produces a 3D plot of the wind velocity magnitude in a specified domain. This returns a Plotly figure for use of downstream visualization.
176+
"""Produces a 3D plot of the wind velocity magnitude in a specified domain. This returns a Plotly figure for use of
177+
downstream visualization.
171178
172179
Parameters
173180
----------
174181
spacing : tuple[float, float, float]
175-
Spacing array that determines the spacing of points to be used in each dimension of the 3D field. Typically, of the form grid_dimensions (a 3x1 vector representing the dimensions of the domain) divided by the grid_levels, which determine the resolution of the wind field in each respective dimension.
182+
Spacing array that determines the spacing of points to be used in each dimension of the 3D field. Typically, of
183+
the form grid_dimensions (a 3x1 vector representing the dimensions of the domain) divided by the grid_levels,
184+
which determine the resolution of the wind field in each respective dimension.
176185
wind_field : np.ndarray
177-
3D wind field, typically of shape (Nx, Ny, Nz, 3) (not C-layout, to be reshaped).
186+
3D wind field, typically of shape (Nx, Ny, Nz, 3) (not C-layout, to be reshaped).
178187
surf_count : int, optional
179188
Number of surfaces to be used, by default 75
180189
reshape : bool, optional
181-
Whether to re-format the given wind field into C-order, typically the desirable choice to match the order of entries of the wind field and the provided spacing, by default True
190+
Whether to re-format the given wind field into C-order, typically the desirable choice to match the order of
191+
entries of the wind field and the provided spacing, by default True
182192
transparent : bool, optional
183-
Whether to set the background of the plot to a transparent background, which renders the same on different backgrounds on which this ``Figure`` could be embedded.
193+
Whether to set the background of the plot to a transparent background, which renders the same on different
194+
backgrounds on which this ``Figure`` could be embedded.
184195
185196
Returns
186197
-------
@@ -192,11 +203,7 @@ def plot_velocity_magnitude(
192203

193204
formatted_wind_field = format_wind_field(wind_field) if reshape else wind_field
194205

195-
wind_magnitude = np.sqrt(
196-
formatted_wind_field[0] ** 2
197-
+ formatted_wind_field[1] ** 2
198-
+ formatted_wind_field[2] ** 2
199-
)
206+
wind_magnitude = np.sqrt(formatted_wind_field[0] ** 2 + formatted_wind_field[1] ** 2 + formatted_wind_field[2] ** 2)
200207

201208
fig = go.Figure(
202209
data=go.Volume(

0 commit comments

Comments
 (0)