From 30e214a87103442946bf9374f47ab0a9d27a22bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Schm=C3=B6lder?= Date: Mon, 2 Dec 2024 10:33:30 +0100 Subject: [PATCH 01/13] [WIP] Add field --- CADETPythonSimulator/field.py | 950 ++++++++++++++++++++++++++++++++++ 1 file changed, 950 insertions(+) create mode 100644 CADETPythonSimulator/field.py diff --git a/CADETPythonSimulator/field.py b/CADETPythonSimulator/field.py new file mode 100644 index 0000000..4d2cfef --- /dev/null +++ b/CADETPythonSimulator/field.py @@ -0,0 +1,950 @@ +from typing import NoReturn, Optional, Union + +import numpy as np +import numpy.typing as npt +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from scipy.interpolate import RegularGridInterpolator + + +class Field: + """ + Represents a single field in a state block. + + Parameters + ---------- + name : str + Name of the field (e.g., "concentration"). + dimensions : dict[str, npt.ArrayLike] + Mapping of dimension names to their coordinate arrays. + n_components : Optional[int], optional + Number of components in the field (e.g., species in a concentration field). + If None, the field is treated as a scalar field. + data : Optional[npt.ArrayLike], optional + Initial data for the field. If None, data is initialized as zeros. + """ + + def __init__( + + self, + name: str, + dimensions: dict[str, npt.ArrayLike], + n_components: Optional[int] = None, + data: Optional[npt.ArrayLike] = None, + ): + self.name = name + self.dimensions = {k: np.array(v) for k, v in dimensions.items()} + self.n_components = n_components + + # Initialize the data array + if data is None: + data = np.zeros(self.shape) + self.data = data + + @property + def is_scalar_field(self) -> bool: + """bool: True if field does not have components, False otherwise.""" + return self.n_components is None + + @property + def data(self) -> np.ndarray: + """np.ndarray: The data array for the field.""" + return self._data + + @data.setter + def data(self, value: npt.ArrayLike): + """Set the data array after validating its shape.""" + value = np.array(value) + if value.shape != self.shape: + raise ValueError( + f"Assigned data shape {value.shape} does not match the expected shape " + f"{self.shape}." + ) + self._data = value + + @property + def n_dimensions(self) -> int: + """int: Return the number of dimensions.""" + return len(self.dimensions) + + @property + def n_cells(self) -> int: + """int: Return the total number of cells from the product of dimensions.""" + return int(np.prod(self.dimension_shape)) + + @property + def dimension_shape(self) -> tuple[int, ...]: + """tuple[int]: Return the shape derived from dimensions.""" + return tuple(len(v) for v in self.dimensions.values()) + + @property + def shape(self) -> tuple[int, ...]: + """tuple[int]: Return the complete shape of the data array.""" + shape = self.dimension_shape + if not self.is_scalar_field: + shape += (self.n_components,) + return shape + + @property + def n_dof(self) -> int: + """int: Return the total number of degrees of freedom.""" + return np.prod(self.shape) + + @property + def data_flat(self) -> np.ndarray: + """np.ndarray: Return the data array flattened into one dimension.""" + return self.data.reshape(-1) + + @data_flat.setter + def data_flat(self, data_flat: np.ndarray) -> None: + """Set the data array from a flattened array.""" + if data_flat.size != self.n_dof: + raise ValueError( + f"Flattened data size {data_flat.size} does not match the " + f"field's degrees of freedom {self.n_dof}." + ) + self.data = data_flat.reshape(self.shape) + + def plot( + self, + title: str = None, + fixed_dims: Optional[dict[str, Union[float, int]]] = None, + method: str = "surface", + fig: Optional[plt.Figure] = None, + axes: Optional[plt.Axes | list[plt.Axes]] = None, + ) -> tuple[plt.Figure, plt.Axes | list[plt.Axes]]: + """ + Plot the field data, supporting 1D, 2D, and slicing for higher dimensions. + + Parameters + ---------- + title : str, optional + Title for the plot. + fixed_dims : dict[str, Union[float, int]], optional + Fixed dimensions to slice for higher-dimensional fields. + Only required for fields with more than 2 dimensions. + method : str, optional + Plotting method: "surface" or "heatmap". Default is "surface". + fig : matplotlib.figure.Figure, optional + Predefined figure to draw on. If not provided, a new figure will be created. + axes : Optional[plt.Axes | list[plt.Axes]], optional + Predefined axes to draw on. If not provided, new axes will be created. + + Returns + ------- + tuple[plt.Figure, plt.Axes | list[plt.Axes]] + The figure and list of axes used for the plot. + + Raises + ------ + ValueError + If the field has more than 2 dimensions and `fixed_dims` is not provided. + If an unknown plotting method is specified. + """ + n_dims = len(self.dimensions) + + # Handle higher-dimensional fields + if n_dims > 2: + if fixed_dims is None: + raise ValueError( + f"Field has {n_dims} dimensions. " + "Please provide fixed_dims to reduce to 1D or 2D for plotting." + ) + # Slice along specified dimensions + sliced_field = self[fixed_dims] + # Recursively call plot on the reduced field + return sliced_field.plot(title=title, method=method, fig=fig, axes=axes) + + # 1D Plot + if n_dims == 1: + x = list(self.dimensions.values())[0] + y = self.data + component_names = ( + [f"Component {i}" for i in range(self.n_components)] + if self.n_components else ["Value"] + ) + return plot_1D(x, y, component_names, title or self.name, fig, axes) + + # 2D Plot + if n_dims == 2: + x, y = self.dimensions.values() + z = self.data + component_names = ( + [f"Component {i}" for i in range(self.n_components)] + if self.n_components else ["Value"] + ) + if method == "surface": + return plot_2D_surface( + x, y, z, component_names, title or self.name, fig, axes + ) + elif method == "heatmap": + return plot_2D_heatmap( + x, y, z, component_names, title or self.name, fig, axes + ) + else: + raise ValueError(f"Unknown plotting method: {method}") + + def resample( + self, + resample_dims: dict[str, Union[int, npt.ArrayLike]] + ) -> "Field": + """ + Resample the field data to new dimensions using interpolation. + + Parameters + ---------- + resample_dims : dict[str, Union[int, npt.ArrayLike]] + Dictionary specifying the resampling for each dimension. + - If a dimension key is missing, it remains unchanged. + - If a value is an integer, it generates a linspace with that many points. + - If a value is an array, it is used directly. + + Returns + ------- + Field + A new Field object with the interpolated data and updated dimensions. + """ + # Validate resample_dims keys + invalid_keys = set(resample_dims.keys()) - set(self.dimensions.keys()) + if invalid_keys: + raise ValueError( + f"Invalid dimension keys in resample_dims: {invalid_keys}. " + f"Valid dimensions are {list(self.dimensions.keys())}." + ) + + # Generate new dimensions + new_dimensions = {} + for dim, original_coords in self.dimensions.items(): + if dim in resample_dims: + value = resample_dims[dim] + if isinstance(value, int): + new_dimensions[dim] = np.linspace( + original_coords.min(), original_coords.max(), value + ) + elif isinstance(value, np.ndarray): + new_dimensions[dim] = value + else: + raise TypeError( + f"Dimension '{dim}' must be an integer or an array." + ) + else: + # Keep the original coordinates if not specified in resample_dims + new_dimensions[dim] = original_coords + + # Use FieldInterpolator for interpolation + interpolated_field = FieldInterpolator(self) + interpolated_data = interpolated_field(**new_dimensions) + + # Return a new Field + return Field( + name=f"{self.name}_resampled", + dimensions=new_dimensions, + n_components=self.n_components, + data=interpolated_data, + ) + + def normalize(self) -> "NormalizedField": + """ + Normalize the field data. + + Returns + ------- + NormalizedField + A new NormalizedField instance with normalized data. + + Notes + ----- + The normalization is performed using the formula: + normalized_data = (data - min(data)) / (max(data) - min(data)) + """ + return NormalizedField(self) + + def derivative(self, dimension: str) -> "Field": + """ + Compute the derivative of the field along a specified dimension. + + Parameters + ---------- + dimension : str + The dimension along which to compute the derivative. + + Returns + ------- + Field + A new Field representing the derivative. + """ + if dimension not in self.dimensions: + raise ValueError(f"Dimension '{dimension}' not found in field.") + + # Use FieldInterpolator for derivative computation + interpolated_field = FieldInterpolator(self) + return interpolated_field.derivative(dimension) + + def anti_derivative(self, dimension: str, initial_value: float = 0) -> "Field": + """ + Compute the anti-derivative of the field along a specified dimension. + + Parameters + ---------- + dimension : str + The dimension along which to compute the anti-derivative. + initial_value : float, optional + Initial value for the anti-derivative. Default is 0. + + Returns + ------- + Field + A new Field representing the anti-derivative. + """ + if dimension not in self.dimensions: + raise ValueError(f"Dimension '{dimension}' not found in field.") + + # Use FieldInterpolator for anti-derivative computation + interpolated_field = FieldInterpolator(self) + return interpolated_field.anti_derivative(dimension, initial_value) + + def __getitem__( + self, + slices: dict[str, float | int] + ) -> Union["Field", np.ndarray]: + """ + Slice the field along specified dimensions. + + Parameters + ---------- + slices : dict[str, float | int] + A dictionary specifying the dimension(s) to slice. + + Returns + ------- + Union[Field, np.ndarray] + A new Field object with reduced dimensions if some dimensions are unspecified. + A scalar or vector (np.ndarray) if all dimensions are specified. + """ + # Ensure valid dimensions are being sliced + for dim in slices.keys(): + if dim not in self.dimensions: + raise KeyError(f"Invalid dimension '{dim}'.") + + # Generate slicing indices + slice_indices = [] + reduced_dimensions = {} + for dim, coords in self.dimensions.items(): + if dim in slices: + # Find the closest index for the specified coordinate + coord = slices[dim] + index = np.searchsorted(coords, coord) + if not (0 <= index < len(coords)): + raise ValueError( + f"Coordinate {coord} out of bounds for dimension '{dim}'." + ) + slice_indices.append(index) + else: + # Keep the dimension if not being sliced + slice_indices.append(slice(None)) + reduced_dimensions[dim] = coords + + # Slice the data + sliced_data = self.data[tuple(slice_indices)] + + # Always return a Field, even when all dimensions are sliced + return Field( + name=self.name, + dimensions=reduced_dimensions, + n_components=self.n_components, + data=sliced_data, + ) + def __repr__(self) -> str: + components = ( + f", components={self.n_components}" + if self.n_components else ", scalar=True" + ) + return f"Field(name='{self.name}', dimensions={list(self.dimensions.keys())}{components})" + + +class NormalizedField(Field): + """ + Represents a normalized version of a Field. + + Parameters + ---------- + field : Field + The original field instance to normalize. + """ + + def __init__(self, field: Field): + # Initialize the base Field class with the same parameters + super().__init__( + name=f"Normalized_{field.name}", + dimensions=field.dimensions, + n_components=field.n_components, + data=field.data, + ) + self.field = field # Keep reference to the original field + self._compute_normalization() + + def _compute_normalization(self): + """Compute the min and max values for normalization.""" + data_min = np.min(self.field.data) + data_max = np.max(self.field.data) + + if data_max == data_min: + raise ValueError("Cannot normalize a field with zero data range.") + + self.data = (self.field.data - data_min) / (data_max - data_min) + + +class FieldInterpolator: + """ + Wrapper around a Field that provides interpolation functionality. + + Parameters + ---------- + field : Field + The original field instance to wrap. + """ + + def __init__(self, field): + self.field = field + self.data = np.array(field.data) # Immutable snapshot of field data + self._initialize_interpolators() + + def _initialize_interpolators(self): + """Initialize interpolators based on the wrapped field's data.""" + grid_points = list(self.field.dimensions.values()) + if self.field.n_components: + self._interpolators = [ + RegularGridInterpolator( + grid_points, self.field.data[..., i], method='pchip' + ) + for i in range(self.field.n_components) + ] + else: + self._interpolator = RegularGridInterpolator(grid_points, self.field.data) + + def __call__(self, **kwargs: npt.ArrayLike) -> np.ndarray: + """ + Evaluate the interpolated field at given coordinates. + + Parameters + ---------- + **kwargs : npt.ArrayLike + Coordinate arrays for each dimension. If a dimension is not specified, + its full original structure will be retained. + + Returns + ------- + np.ndarray + Interpolated values at the specified coordinates. + If some dimensions are unspecified, their original structure is preserved. + """ + # Separate provided and missing dimensions + provided_dims = set(kwargs.keys()) + all_dims = set(self.field.dimensions.keys()) + missing_dims = all_dims - provided_dims + + # Validate dimensions + if not provided_dims.issubset(all_dims): + invalid_dims = provided_dims - all_dims + raise ValueError(f"Invalid dimensions specified: {invalid_dims}") + + # Create query points, filling in original coordinates for missing dimensions + query_coords = [ + kwargs.get(dim, self.field.dimensions[dim]) + for dim in self.field.dimensions.keys() + ] + mesh = np.meshgrid(*query_coords, indexing="ij") + query_points = np.stack([grid.ravel() for grid in mesh], axis=-1) + + # Perform interpolation + if self.field.n_components: + interpolated_components = [ + interp(query_points).reshape(mesh[0].shape) + for interp in self._interpolators + ] + result = np.stack(interpolated_components, axis=-1) + else: + result = self._interpolator(query_points).reshape(mesh[0].shape) + + # Reshape the result based on the query to retain missing dimensions + output_shape = self._determine_output_shape(kwargs) + return result.reshape(output_shape) + + def _determine_output_shape( + self, + query: dict[str, Union[float, npt.ArrayLike]] + ) -> tuple[int, ...]: + """ + Determine the shape of the interpolated result based on the query. + + Parameters + ---------- + query : dict[str, Union[float, npt.ArrayLike]] + The coordinates provided for each dimension in the query. + + Returns + ------- + tuple[int, ...] + The shape of the interpolated result. + """ + output_shape = [] + for dim, original_coords in self.field.dimensions.items(): + if dim in query: + queried_value = query[dim] + if np.isscalar(queried_value): + continue # Scalars do not contribute to the shape + else: + output_shape.append(len(queried_value)) # Retain array length + else: + output_shape.append(len(original_coords)) # Retain full original dimension + + # Add n_components as the last dimension for vector fields + if self.field.n_components: + output_shape.append(self.field.n_components) + + return tuple(output_shape) + + def derivative(self, dimension: str): + """ + Compute the derivative of the field along the specified dimension. + + Parameters + ---------- + dimension : str + Dimension along which to compute the derivative. + + Returns + ------- + Field + A new field representing the derivative. + """ + if dimension not in self.field.dimensions: + raise ValueError(f"Dimension '{dimension}' not found in field.") + + axis = list(self.field.dimensions.keys()).index(dimension) + coords = self.field.dimensions[dimension] + + # Compute derivative for multi-component fields + if self.field.n_components: + derivative_data = np.stack([ + np.gradient(self.field.data[..., i], coords, axis=axis) + for i in range(self.field.n_components) + ], axis=-1) + else: + derivative_data = np.gradient(self.field.data, coords, axis=axis) + + # Return new Field for the derivative + return Field( + name=f"{self.field.name}_derivative_{dimension}", + dimensions=self.field.dimensions, + n_components=self.field.n_components, + data=derivative_data + ) + + def anti_derivative(self, dimension: str, initial_value: float = 0): + """ + Compute the anti-derivative of the field along the specified dimension. + + Parameters + ---------- + dimension : str + Dimension along which to compute the anti-derivative. + initial_value : float, optional + Initial value for the anti-derivative. Default is 0. + + Returns + ------- + Field + A new field representing the anti-derivative. + """ + if dimension not in self.field.dimensions: + raise ValueError(f"Dimension '{dimension}' not found in field.") + + axis = list(self.field.dimensions.keys()).index(dimension) + coords = self.field.dimensions[dimension] + + # Compute anti-derivative for multi-component fields + if self.field.n_components: + anti_derivative_data = np.stack([ + np.cumsum( + self.field.data[..., i] * np.gradient(coords), axis=axis + ) + initial_value + for i in range(self.field.n_components) + ], axis=-1) + else: + anti_derivative_data = np.cumsum( + self.field.data * np.gradient(coords), axis=axis + ) + initial_value + + # Return new Field for the anti-derivative + return Field( + name=f"{self.field.name}_anti_derivative_{dimension}", + dimensions=self.field.dimensions, + n_components=self.field.n_components, + data=anti_derivative_data + ) + + def __repr__(self) -> str: + return f"FieldInterpolator({self.field})" + + +def plot_1D( + x: np.ndarray, + y: np.ndarray, + component_names: list[str], + title: str, + fig: Optional[plt.Figure] = None, + ax: Optional[plt.Axes] = None, + ) -> tuple[plt.Figure, plt.Axes]: + """ + Plot a 1D field with a line for each component. + + Parameters + ---------- + x : np.ndarray + The x-axis data (1D coordinate array). + y : np.ndarray + The y-axis data, with shape (len(x), n_components). + component_names : list[str] + Names of the components to label in the legend. + title : str + Title of the plot. + fig : matplotlib.figure.Figure, optional + Predefined figure to draw on. If not provided, a new figure will be created. + ax : matplotlib.axes.Axes, optional + Predefined axis to draw on. If not provided, a new axis will be created. + + Returns + ------- + tuple[matplotlib.figure.Figure, list[matplotlib.axes.Axes]] + The figure and list containing the single axis used for the plot. + """ + if fig is None or ax is None: + fig, ax = plt.subplots() + + for i, component in enumerate(component_names): + ax.plot(x, y[:, i], label=component) + + ax.set_title(title) + ax.set_xlabel("X-axis") + ax.set_ylabel("Values") + ax.legend() + + return fig, ax + + +def plot_2D_surface( + x: np.ndarray, + y: np.ndarray, + z: np.ndarray, + component_names: list[str], + title: str, + fig: Optional[plt.Figure] = None, + axes: Optional[list[plt.Axes]] = None, + ) -> tuple[plt.Figure, list[plt.Axes]]: + """ + Plot a 2D field with a surface for each component. + + Parameters + ---------- + x : np.ndarray + The x-axis data (1D coordinate array). + y : np.ndarray + The y-axis data (1D coordinate array). + z : np.ndarray + The z-axis data, with shape (len(x), len(y), n_components). + component_names : list[str] + Names of the components to label the subplots. + title : str + Title of the plot. + fig : matplotlib.figure.Figure, optional + Predefined figure to draw on. If not provided, a new figure will be created. + axes : list[matplotlib.axes.Axes], optional + Predefined axes to draw on. If not provided, new axes will be created. + + Returns + ------- + tuple[matplotlib.figure.Figure, list[matplotlib.axes.Axes]] + The figure and list of axes used for the plot. + """ + n_components = z.shape[-1] + + if axes is None: + fig, axes = plt.subplots( + 1, n_components, + figsize=(5 * n_components, 5), + subplot_kw={'projection': '3d'} + ) + axes = [axes] if n_components == 1 else axes + + X, Y = np.meshgrid(x, y, indexing="ij") + for i, ax in enumerate(axes): + ax.plot_surface(X, Y, z[..., i], cmap="viridis") + ax.set_title(f"{title} - {component_names[i]}") + ax.set_xlabel("X-axis") + ax.set_ylabel("Y-axis") + ax.set_zlabel("Values") + + return fig, axes + + +def plot_2D_heatmap( + x: np.ndarray, + y: np.ndarray, + z: np.ndarray, + component_names: list[str], + title: str, + fig: Optional[plt.Figure] = None, + axes: Optional[list[plt.Axes]] = None, + ) -> tuple[plt.Figure, list[plt.Axes]]: + """ + Plot a 2D field with a heatmap for each component. + + Parameters + ---------- + x : np.ndarray + The x-axis data (1D coordinate array). + y : np.ndarray + The y-axis data (1D coordinate array). + z : np.ndarray + The z-axis data, with shape (len(x), len(y), n_components). + component_names : list[str] + Names of the components to label the subplots. + title : str + Title of the plot. + fig : matplotlib.figure.Figure, optional + Predefined figure to draw on. If not provided, a new figure will be created. + axes : list[matplotlib.axes.Axes], optional + Predefined axes to draw on. If not provided, new axes will be created. + + Returns + ------- + tuple[matplotlib.figure.Figure, list[matplotlib.axes.Axes]] + The figure and list of axes used for the plot. + """ + n_components = z.shape[-1] + + if axes is None: + fig, axes = plt.subplots( + 1, n_components, + figsize=(5 * n_components, 5), + constrained_layout=True + ) + axes = [axes] if n_components == 1 else axes + + X, Y = np.meshgrid(x, y, indexing="ij") + for i, ax in enumerate(axes): + c = ax.pcolormesh(X, Y, z[..., i], cmap="viridis", shading="auto") + fig.colorbar(c, ax=ax, orientation="vertical") + ax.set_title(f"{title} - {component_names[i]}") + ax.set_xlabel("X-axis") + ax.set_ylabel("Y-axis") + + return fig, axes + + +# %% Testing utilities + +def assert_equal(value, expected, message=""): + message = f"Test failed: {message}. Expected {expected}, got {value}." + assert value == expected, message + + +def assert_shape(shape, expected_shape, message=""): + message = f"Test failed: {message}. Expected shape {expected_shape}, got {shape}." + assert shape == expected_shape, message + + +# %% Initialization + +def test_field_initialization(): + # Scalar field + dimensions = { + "axial": np.linspace(0, 10, 11), + "radial": np.linspace(0, 5, 6) + } + viscosity = Field(name="viscosity", dimensions=dimensions) + assert_shape(viscosity.shape, (11, 6), "Scalar field shape") + assert_equal(viscosity.n_dof, 11 * 6, "Scalar field degrees of freedom") + + # Vector field + concentration = Field(name="concentration", dimensions=dimensions, n_components=3) + assert_shape(concentration.shape, (11, 6, 3), "Vector field shape") + assert_equal(concentration.n_dof, 11 * 6 * 3, "Vector field degrees of freedom") + + # Custom data + data = np.ones((11, 6, 3)) + concentration_with_data = Field(name="concentration", dimensions=dimensions, n_components=3, data=data) + assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") + + +# %% Plotting + +def test_plotting(): + # 1D Plot + dimensions = {"x": np.linspace(0, 10, 11)} + field_1D = Field(name="1D Field", dimensions=dimensions, n_components=2, data=np.random.random((11, 2))) + fig, ax = field_1D.plot() + assert isinstance(ax, plt.Axes), "1D plot returns one axis" + + # 2D Plot + dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} + field_2D = Field(name="2D Field", dimensions=dimensions, n_components=3, data=np.random.random((11, 6, 3))) + fig, axes = field_2D.plot() + assert len(axes) == 3, "2D plot returns one axis per component" + + +# %% Slicing + +def test_field_slicing(): + dimensions = {"axial": np.linspace(0, 10, 11), "radial": np.linspace(0, 5, 6)} + field = Field(name="concentration", dimensions=dimensions, n_components=3) + + # Slice along one dimension + field_sliced = field[{"axial": 0}] + assert_equal(len(field_sliced.dimensions), 1, "Field slicing reduces dimensionality") + assert_shape(field_sliced.shape, (6, 3), "Field slicing shape") + + # Slice along all dimensions + field_sliced_all = field[{"axial": 0, "radial": 0}] + assert_equal(len(field_sliced_all.dimensions), 0, "Full slicing removes all dimensions") + assert_shape(field_sliced_all.shape, (3,), "Full slicing results in vector") + + +# %% Normalization + +def test_field_normalization(): + # Define dimensions and data + x = np.linspace(0, 10, 11) + y = np.linspace(0, 5, 6) + z = np.outer(np.sin(x), np.cos(y)) # Example data + + # Create a field + field = Field(name="temperature", dimensions={"x": x, "y": y}, data=z) + + # Normalize the field + normalized_field = field.normalize() + + # Test 1: Check dimensions and structure + assert field.shape == normalized_field.shape, "Normalized field shape mismatch." + np.testing.assert_equal(field.dimensions, normalized_field.dimensions) + + # Test 2: Verify data normalization + normalized_data = normalized_field.data + assert np.isclose(np.min(normalized_data), 0.0), "Normalized data minimum is not 0." + assert np.isclose(np.max(normalized_data), 1.0), "Normalized data maximum is not 1." + + # Test 3: Ensure original field is unchanged + assert np.array_equal(field.data, z), "Original field data was modified." + + + +# %% Interpolation and Resampling + +def test_temperature_use_case(): + temperature_profile_dimensions = { + "time": np.linspace(0, 3600, 3601), + "axial": np.linspace(0, 0.5, 6), + } + + # Idea, only change over time, start with T = 20C, end with 30C + temperature_data = 20 * np.ones((3601, 6)) + + temperature_field = Field( + name="temperature", + dimensions=temperature_profile_dimensions, + data=temperature_data, + ) + + temperature_interpolator = FieldInterpolator(temperature_field) + + def calculate_temperature_at_t_x(t, x): + return temperature_interpolator(time=t, axial=x) + + def calculate_adsorption_from_temperature(k_0, k_1, T): + return k_0 * np.exp(k_1 / T) + + def calculate_parameter(t, x): + T = calculate_temperature_at_t_x(t, x) + k = calculate_adsorption_from_temperature(k_0, k_1, T) + + def model_equation(t): + ... + + +def test_interpolated_field(): + dimensions = { + "axial": np.linspace(0, 10, 11), + "radial": np.linspace(0, 5, 6) + } + data = np.random.random((11, 6, 3)) + concentration = Field(name="concentration", dimensions=dimensions, n_components=3, data=data) + + # Interpolated field + interp_field = FieldInterpolator(concentration) + result = interp_field(axial=1.5, radial=2.1) + assert_shape(result.shape, (3,), "FieldInterpolator output components") + + +def test_resampling(): + dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} + field = Field(name="concentration", dimensions=dimensions, n_components=2, data=np.random.random((11, 6, 2))) + + # Resample one dimension + resampled_field = field.resample({"x": 50}) + assert_shape(resampled_field.shape, (50, 6, 2), "Resampling one dimension") + + # Resample all dimensions + resampled_field_all = field.resample({"x": 50, "y": 25}) + assert_shape(resampled_field_all.shape, (50, 25, 2), "Resampling all dimensions") + + +def test_field_interpolation_and_derivatives(): + # Define test dimensions and data + x = np.linspace(0, 10, 11) # 11 points along x + y = np.linspace(0, 5, 6) # 6 points along y + z = np.outer(np.sin(x), np.cos(y)) # Generate 2D scalar field data + + # Create a Field + field = Field(name="test_field", dimensions={"x": x, "y": y}, data=z) + + # Wrap with FieldInterpolator + interp_field = FieldInterpolator(field) + + # Test 1: Evaluate at an arbitrary point + eval_result = interp_field(x=2.5) + assert_shape( + eval_result.shape, (6, ), + "Evaluation result should be a scalar for scalar fields." + ) + eval_result = interp_field(x=2.5, y=1.1) + assert_shape( + eval_result.shape, (), + "Evaluation result should be a scalar for scalar fields." + ) + + # Test 2: Compute derivative along x + dx_field = interp_field.derivative("x") + assert dx_field.data.shape == field.data.shape, "Derivative shape mismatch!" + + # Test 3: Compute anti-derivative along y + int_y_field = interp_field.anti_derivative("y", initial_value=0) + assert int_y_field.data.shape == field.data.shape, "Anti-derivative shape mismatch!" + + # Test 4: Verify dimensions + assert dx_field.data.shape == field.data.shape, "Derivative shape mismatch!" + assert int_y_field.data.shape == field.data.shape, "Anti-derivative shape mismatch!" + + # Test 5: Edge Case - Evaluate at boundary + boundary_value = interp_field(x=0, y=5) + assert_shape( + boundary_value.shape, (), + "Boundary evaluation should return a scalar for scalar fields." + ) + + +# %% Run tests + +import pytest +if __name__ == "__main__": + pytest.main([str(__file__)]) From 213c3d6e45213acd62bea1d5636518ad3cdc6edb Mon Sep 17 00:00:00 2001 From: daklauss Date: Mon, 2 Dec 2024 11:11:36 +0100 Subject: [PATCH 02/13] Fix some descriptions for ruff --- CADETPythonSimulator/field.py | 32 +++++++++++++++++++++++++++----- 1 file changed, 27 insertions(+), 5 deletions(-) diff --git a/CADETPythonSimulator/field.py b/CADETPythonSimulator/field.py index 4d2cfef..478c473 100644 --- a/CADETPythonSimulator/field.py +++ b/CADETPythonSimulator/field.py @@ -22,16 +22,17 @@ class Field: If None, the field is treated as a scalar field. data : Optional[npt.ArrayLike], optional Initial data for the field. If None, data is initialized as zeros. + """ def __init__( - self, name: str, dimensions: dict[str, npt.ArrayLike], n_components: Optional[int] = None, data: Optional[npt.ArrayLike] = None, ): + """Construct the object.""" self.name = name self.dimensions = {k: np.array(v) for k, v in dimensions.items()} self.n_components = n_components @@ -140,6 +141,7 @@ def plot( ValueError If the field has more than 2 dimensions and `fixed_dims` is not provided. If an unknown plotting method is specified. + """ n_dims = len(self.dimensions) @@ -203,6 +205,7 @@ def resample( ------- Field A new Field object with the interpolated data and updated dimensions. + """ # Validate resample_dims keys invalid_keys = set(resample_dims.keys()) - set(self.dimensions.keys()) @@ -256,6 +259,7 @@ def normalize(self) -> "NormalizedField": ----- The normalization is performed using the formula: normalized_data = (data - min(data)) / (max(data) - min(data)) + """ return NormalizedField(self) @@ -272,6 +276,7 @@ def derivative(self, dimension: str) -> "Field": ------- Field A new Field representing the derivative. + """ if dimension not in self.dimensions: raise ValueError(f"Dimension '{dimension}' not found in field.") @@ -295,6 +300,7 @@ def anti_derivative(self, dimension: str, initial_value: float = 0) -> "Field": ------- Field A new Field representing the anti-derivative. + """ if dimension not in self.dimensions: raise ValueError(f"Dimension '{dimension}' not found in field.") @@ -318,8 +324,10 @@ def __getitem__( Returns ------- Union[Field, np.ndarray] - A new Field object with reduced dimensions if some dimensions are unspecified. + A new Field object with reduced dimensions if some dimensions are + unspecified. A scalar or vector (np.ndarray) if all dimensions are specified. + """ # Ensure valid dimensions are being sliced for dim in slices.keys(): @@ -355,11 +363,13 @@ def __getitem__( data=sliced_data, ) def __repr__(self) -> str: + """Represantator.""" components = ( f", components={self.n_components}" if self.n_components else ", scalar=True" ) - return f"Field(name='{self.name}', dimensions={list(self.dimensions.keys())}{components})" + return f"Field(name='{self.name}')"\ + +f", dimensions={list(self.dimensions.keys())}{components}" class NormalizedField(Field): @@ -370,10 +380,11 @@ class NormalizedField(Field): ---------- field : Field The original field instance to normalize. + """ def __init__(self, field: Field): - # Initialize the base Field class with the same parameters + """Initialize the base Field class with the same parameters.""" super().__init__( name=f"Normalized_{field.name}", dimensions=field.dimensions, @@ -402,9 +413,11 @@ class FieldInterpolator: ---------- field : Field The original field instance to wrap. + """ def __init__(self, field): + """Construct the object.""" self.field = field self.data = np.array(field.data) # Immutable snapshot of field data self._initialize_interpolators() @@ -437,6 +450,7 @@ def __call__(self, **kwargs: npt.ArrayLike) -> np.ndarray: np.ndarray Interpolated values at the specified coordinates. If some dimensions are unspecified, their original structure is preserved. + """ # Separate provided and missing dimensions provided_dims = set(kwargs.keys()) @@ -486,6 +500,7 @@ def _determine_output_shape( ------- tuple[int, ...] The shape of the interpolated result. + """ output_shape = [] for dim, original_coords in self.field.dimensions.items(): @@ -496,7 +511,8 @@ def _determine_output_shape( else: output_shape.append(len(queried_value)) # Retain array length else: - output_shape.append(len(original_coords)) # Retain full original dimension + output_shape.append(len(original_coords)) + # Retain full original dimension # Add n_components as the last dimension for vector fields if self.field.n_components: @@ -517,6 +533,7 @@ def derivative(self, dimension: str): ------- Field A new field representing the derivative. + """ if dimension not in self.field.dimensions: raise ValueError(f"Dimension '{dimension}' not found in field.") @@ -556,6 +573,7 @@ def anti_derivative(self, dimension: str, initial_value: float = 0): ------- Field A new field representing the anti-derivative. + """ if dimension not in self.field.dimensions: raise ValueError(f"Dimension '{dimension}' not found in field.") @@ -585,6 +603,7 @@ def anti_derivative(self, dimension: str, initial_value: float = 0): ) def __repr__(self) -> str: + """Return represantation.""" return f"FieldInterpolator({self.field})" @@ -618,6 +637,7 @@ def plot_1D( ------- tuple[matplotlib.figure.Figure, list[matplotlib.axes.Axes]] The figure and list containing the single axis used for the plot. + """ if fig is None or ax is None: fig, ax = plt.subplots() @@ -666,6 +686,7 @@ def plot_2D_surface( ------- tuple[matplotlib.figure.Figure, list[matplotlib.axes.Axes]] The figure and list of axes used for the plot. + """ n_components = z.shape[-1] @@ -721,6 +742,7 @@ def plot_2D_heatmap( ------- tuple[matplotlib.figure.Figure, list[matplotlib.axes.Axes]] The figure and list of axes used for the plot. + """ n_components = z.shape[-1] From 481d5f0b18b1fb7c0aabd5107ab81865ce06d3be Mon Sep 17 00:00:00 2001 From: daklauss Date: Mon, 2 Dec 2024 11:14:58 +0100 Subject: [PATCH 03/13] Created test field and moved tests to new file --- CADETPythonSimulator/field.py | 206 --------------------------------- tests/test_field.py | 208 ++++++++++++++++++++++++++++++++++ 2 files changed, 208 insertions(+), 206 deletions(-) create mode 100644 tests/test_field.py diff --git a/CADETPythonSimulator/field.py b/CADETPythonSimulator/field.py index 478c473..8cc784b 100644 --- a/CADETPythonSimulator/field.py +++ b/CADETPythonSimulator/field.py @@ -764,209 +764,3 @@ def plot_2D_heatmap( return fig, axes - -# %% Testing utilities - -def assert_equal(value, expected, message=""): - message = f"Test failed: {message}. Expected {expected}, got {value}." - assert value == expected, message - - -def assert_shape(shape, expected_shape, message=""): - message = f"Test failed: {message}. Expected shape {expected_shape}, got {shape}." - assert shape == expected_shape, message - - -# %% Initialization - -def test_field_initialization(): - # Scalar field - dimensions = { - "axial": np.linspace(0, 10, 11), - "radial": np.linspace(0, 5, 6) - } - viscosity = Field(name="viscosity", dimensions=dimensions) - assert_shape(viscosity.shape, (11, 6), "Scalar field shape") - assert_equal(viscosity.n_dof, 11 * 6, "Scalar field degrees of freedom") - - # Vector field - concentration = Field(name="concentration", dimensions=dimensions, n_components=3) - assert_shape(concentration.shape, (11, 6, 3), "Vector field shape") - assert_equal(concentration.n_dof, 11 * 6 * 3, "Vector field degrees of freedom") - - # Custom data - data = np.ones((11, 6, 3)) - concentration_with_data = Field(name="concentration", dimensions=dimensions, n_components=3, data=data) - assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") - - -# %% Plotting - -def test_plotting(): - # 1D Plot - dimensions = {"x": np.linspace(0, 10, 11)} - field_1D = Field(name="1D Field", dimensions=dimensions, n_components=2, data=np.random.random((11, 2))) - fig, ax = field_1D.plot() - assert isinstance(ax, plt.Axes), "1D plot returns one axis" - - # 2D Plot - dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} - field_2D = Field(name="2D Field", dimensions=dimensions, n_components=3, data=np.random.random((11, 6, 3))) - fig, axes = field_2D.plot() - assert len(axes) == 3, "2D plot returns one axis per component" - - -# %% Slicing - -def test_field_slicing(): - dimensions = {"axial": np.linspace(0, 10, 11), "radial": np.linspace(0, 5, 6)} - field = Field(name="concentration", dimensions=dimensions, n_components=3) - - # Slice along one dimension - field_sliced = field[{"axial": 0}] - assert_equal(len(field_sliced.dimensions), 1, "Field slicing reduces dimensionality") - assert_shape(field_sliced.shape, (6, 3), "Field slicing shape") - - # Slice along all dimensions - field_sliced_all = field[{"axial": 0, "radial": 0}] - assert_equal(len(field_sliced_all.dimensions), 0, "Full slicing removes all dimensions") - assert_shape(field_sliced_all.shape, (3,), "Full slicing results in vector") - - -# %% Normalization - -def test_field_normalization(): - # Define dimensions and data - x = np.linspace(0, 10, 11) - y = np.linspace(0, 5, 6) - z = np.outer(np.sin(x), np.cos(y)) # Example data - - # Create a field - field = Field(name="temperature", dimensions={"x": x, "y": y}, data=z) - - # Normalize the field - normalized_field = field.normalize() - - # Test 1: Check dimensions and structure - assert field.shape == normalized_field.shape, "Normalized field shape mismatch." - np.testing.assert_equal(field.dimensions, normalized_field.dimensions) - - # Test 2: Verify data normalization - normalized_data = normalized_field.data - assert np.isclose(np.min(normalized_data), 0.0), "Normalized data minimum is not 0." - assert np.isclose(np.max(normalized_data), 1.0), "Normalized data maximum is not 1." - - # Test 3: Ensure original field is unchanged - assert np.array_equal(field.data, z), "Original field data was modified." - - - -# %% Interpolation and Resampling - -def test_temperature_use_case(): - temperature_profile_dimensions = { - "time": np.linspace(0, 3600, 3601), - "axial": np.linspace(0, 0.5, 6), - } - - # Idea, only change over time, start with T = 20C, end with 30C - temperature_data = 20 * np.ones((3601, 6)) - - temperature_field = Field( - name="temperature", - dimensions=temperature_profile_dimensions, - data=temperature_data, - ) - - temperature_interpolator = FieldInterpolator(temperature_field) - - def calculate_temperature_at_t_x(t, x): - return temperature_interpolator(time=t, axial=x) - - def calculate_adsorption_from_temperature(k_0, k_1, T): - return k_0 * np.exp(k_1 / T) - - def calculate_parameter(t, x): - T = calculate_temperature_at_t_x(t, x) - k = calculate_adsorption_from_temperature(k_0, k_1, T) - - def model_equation(t): - ... - - -def test_interpolated_field(): - dimensions = { - "axial": np.linspace(0, 10, 11), - "radial": np.linspace(0, 5, 6) - } - data = np.random.random((11, 6, 3)) - concentration = Field(name="concentration", dimensions=dimensions, n_components=3, data=data) - - # Interpolated field - interp_field = FieldInterpolator(concentration) - result = interp_field(axial=1.5, radial=2.1) - assert_shape(result.shape, (3,), "FieldInterpolator output components") - - -def test_resampling(): - dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} - field = Field(name="concentration", dimensions=dimensions, n_components=2, data=np.random.random((11, 6, 2))) - - # Resample one dimension - resampled_field = field.resample({"x": 50}) - assert_shape(resampled_field.shape, (50, 6, 2), "Resampling one dimension") - - # Resample all dimensions - resampled_field_all = field.resample({"x": 50, "y": 25}) - assert_shape(resampled_field_all.shape, (50, 25, 2), "Resampling all dimensions") - - -def test_field_interpolation_and_derivatives(): - # Define test dimensions and data - x = np.linspace(0, 10, 11) # 11 points along x - y = np.linspace(0, 5, 6) # 6 points along y - z = np.outer(np.sin(x), np.cos(y)) # Generate 2D scalar field data - - # Create a Field - field = Field(name="test_field", dimensions={"x": x, "y": y}, data=z) - - # Wrap with FieldInterpolator - interp_field = FieldInterpolator(field) - - # Test 1: Evaluate at an arbitrary point - eval_result = interp_field(x=2.5) - assert_shape( - eval_result.shape, (6, ), - "Evaluation result should be a scalar for scalar fields." - ) - eval_result = interp_field(x=2.5, y=1.1) - assert_shape( - eval_result.shape, (), - "Evaluation result should be a scalar for scalar fields." - ) - - # Test 2: Compute derivative along x - dx_field = interp_field.derivative("x") - assert dx_field.data.shape == field.data.shape, "Derivative shape mismatch!" - - # Test 3: Compute anti-derivative along y - int_y_field = interp_field.anti_derivative("y", initial_value=0) - assert int_y_field.data.shape == field.data.shape, "Anti-derivative shape mismatch!" - - # Test 4: Verify dimensions - assert dx_field.data.shape == field.data.shape, "Derivative shape mismatch!" - assert int_y_field.data.shape == field.data.shape, "Anti-derivative shape mismatch!" - - # Test 5: Edge Case - Evaluate at boundary - boundary_value = interp_field(x=0, y=5) - assert_shape( - boundary_value.shape, (), - "Boundary evaluation should return a scalar for scalar fields." - ) - - -# %% Run tests - -import pytest -if __name__ == "__main__": - pytest.main([str(__file__)]) diff --git a/tests/test_field.py b/tests/test_field.py new file mode 100644 index 0000000..ad6d708 --- /dev/null +++ b/tests/test_field.py @@ -0,0 +1,208 @@ +import numpy as np +from CADETPythonSimulator.field import Field, FieldInterpolator +import matplotlib.pyplot as plt +# %% Testing utilities + +def assert_equal(value, expected, message=""): + message = f"Test failed: {message}. Expected {expected}, got {value}." + assert value == expected, message + + +def assert_shape(shape, expected_shape, message=""): + message = f"Test failed: {message}. Expected shape {expected_shape}, got {shape}." + assert shape == expected_shape, message + + +# %% Initialization + +def test_field_initialization(): + # Scalar field + dimensions = { + "axial": np.linspace(0, 10, 11), + "radial": np.linspace(0, 5, 6) + } + viscosity = Field(name="viscosity", dimensions=dimensions) + assert_shape(viscosity.shape, (11, 6), "Scalar field shape") + assert_equal(viscosity.n_dof, 11 * 6, "Scalar field degrees of freedom") + + # Vector field + concentration = Field(name="concentration", dimensions=dimensions, n_components=3) + assert_shape(concentration.shape, (11, 6, 3), "Vector field shape") + assert_equal(concentration.n_dof, 11 * 6 * 3, "Vector field degrees of freedom") + + # Custom data + data = np.ones((11, 6, 3)) + concentration_with_data = Field(name="concentration", dimensions=dimensions, n_components=3, data=data) + assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") + + +# %% Plotting + +def test_plotting(): + # 1D Plot + dimensions = {"x": np.linspace(0, 10, 11)} + field_1D = Field(name="1D Field", dimensions=dimensions, n_components=2, data=np.random.random((11, 2))) + fig, ax = field_1D.plot() + assert isinstance(ax, plt.Axes), "1D plot returns one axis" + + # 2D Plot + dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} + field_2D = Field(name="2D Field", dimensions=dimensions, n_components=3, data=np.random.random((11, 6, 3))) + fig, axes = field_2D.plot() + assert len(axes) == 3, "2D plot returns one axis per component" + + +# %% Slicing + +def test_field_slicing(): + dimensions = {"axial": np.linspace(0, 10, 11), "radial": np.linspace(0, 5, 6)} + field = Field(name="concentration", dimensions=dimensions, n_components=3) + + # Slice along one dimension + field_sliced = field[{"axial": 0}] + assert_equal(len(field_sliced.dimensions), 1, "Field slicing reduces dimensionality") + assert_shape(field_sliced.shape, (6, 3), "Field slicing shape") + + # Slice along all dimensions + field_sliced_all = field[{"axial": 0, "radial": 0}] + assert_equal(len(field_sliced_all.dimensions), 0, "Full slicing removes all dimensions") + assert_shape(field_sliced_all.shape, (3,), "Full slicing results in vector") + + +# %% Normalization + +def test_field_normalization(): + # Define dimensions and data + x = np.linspace(0, 10, 11) + y = np.linspace(0, 5, 6) + z = np.outer(np.sin(x), np.cos(y)) # Example data + + # Create a field + field = Field(name="temperature", dimensions={"x": x, "y": y}, data=z) + + # Normalize the field + normalized_field = field.normalize() + + # Test 1: Check dimensions and structure + assert field.shape == normalized_field.shape, "Normalized field shape mismatch." + np.testing.assert_equal(field.dimensions, normalized_field.dimensions) + + # Test 2: Verify data normalization + normalized_data = normalized_field.data + assert np.isclose(np.min(normalized_data), 0.0), "Normalized data minimum is not 0." + assert np.isclose(np.max(normalized_data), 1.0), "Normalized data maximum is not 1." + + # Test 3: Ensure original field is unchanged + assert np.array_equal(field.data, z), "Original field data was modified." + + + +# %% Interpolation and Resampling + +def test_temperature_use_case(): + temperature_profile_dimensions = { + "time": np.linspace(0, 3600, 3601), + "axial": np.linspace(0, 0.5, 6), + } + + # Idea, only change over time, start with T = 20C, end with 30C + temperature_data = 20 * np.ones((3601, 6)) + + temperature_field = Field( + name="temperature", + dimensions=temperature_profile_dimensions, + data=temperature_data, + ) + + temperature_interpolator = FieldInterpolator(temperature_field) + + def calculate_temperature_at_t_x(t, x): + return temperature_interpolator(time=t, axial=x) + + def calculate_adsorption_from_temperature(k_0, k_1, T): + return k_0 * np.exp(k_1 / T) + + def calculate_parameter(t, x): + T = calculate_temperature_at_t_x(t, x) + k = calculate_adsorption_from_temperature(k_0, k_1, T) + + def model_equation(t): + ... + + +def test_interpolated_field(): + dimensions = { + "axial": np.linspace(0, 10, 11), + "radial": np.linspace(0, 5, 6) + } + data = np.random.random((11, 6, 3)) + concentration = Field(name="concentration", dimensions=dimensions, n_components=3, data=data) + + # Interpolated field + interp_field = FieldInterpolator(concentration) + result = interp_field(axial=1.5, radial=2.1) + assert_shape(result.shape, (3,), "FieldInterpolator output components") + + +def test_resampling(): + dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} + field = Field(name="concentration", dimensions=dimensions, n_components=2, data=np.random.random((11, 6, 2))) + + # Resample one dimension + resampled_field = field.resample({"x": 50}) + assert_shape(resampled_field.shape, (50, 6, 2), "Resampling one dimension") + + # Resample all dimensions + resampled_field_all = field.resample({"x": 50, "y": 25}) + assert_shape(resampled_field_all.shape, (50, 25, 2), "Resampling all dimensions") + + +def test_field_interpolation_and_derivatives(): + # Define test dimensions and data + x = np.linspace(0, 10, 11) # 11 points along x + y = np.linspace(0, 5, 6) # 6 points along y + z = np.outer(np.sin(x), np.cos(y)) # Generate 2D scalar field data + + # Create a Field + field = Field(name="test_field", dimensions={"x": x, "y": y}, data=z) + + # Wrap with FieldInterpolator + interp_field = FieldInterpolator(field) + + # Test 1: Evaluate at an arbitrary point + eval_result = interp_field(x=2.5) + assert_shape( + eval_result.shape, (6, ), + "Evaluation result should be a scalar for scalar fields." + ) + eval_result = interp_field(x=2.5, y=1.1) + assert_shape( + eval_result.shape, (), + "Evaluation result should be a scalar for scalar fields." + ) + + # Test 2: Compute derivative along x + dx_field = interp_field.derivative("x") + assert dx_field.data.shape == field.data.shape, "Derivative shape mismatch!" + + # Test 3: Compute anti-derivative along y + int_y_field = interp_field.anti_derivative("y", initial_value=0) + assert int_y_field.data.shape == field.data.shape, "Anti-derivative shape mismatch!" + + # Test 4: Verify dimensions + assert dx_field.data.shape == field.data.shape, "Derivative shape mismatch!" + assert int_y_field.data.shape == field.data.shape, "Anti-derivative shape mismatch!" + + # Test 5: Edge Case - Evaluate at boundary + boundary_value = interp_field(x=0, y=5) + assert_shape( + boundary_value.shape, (), + "Boundary evaluation should return a scalar for scalar fields." + ) + + +# %% Run tests + +import pytest +if __name__ == "__main__": + pytest.main('test_field.py') From 08a5f9ddfa7d855b2cd5eb888b42f22772c83229 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Sun, 8 Dec 2024 12:45:53 +0100 Subject: [PATCH 04/13] Improve Field test coverage (#22) --- tests/test_field.py | 52 ++++++++++++++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 17 deletions(-) diff --git a/tests/test_field.py b/tests/test_field.py index ad6d708..717867e 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -1,8 +1,10 @@ +import pytest import numpy as np from CADETPythonSimulator.field import Field, FieldInterpolator import matplotlib.pyplot as plt # %% Testing utilities + def assert_equal(value, expected, message=""): message = f"Test failed: {message}. Expected {expected}, got {value}." assert value == expected, message @@ -26,14 +28,21 @@ def test_field_initialization(): assert_equal(viscosity.n_dof, 11 * 6, "Scalar field degrees of freedom") # Vector field - concentration = Field(name="concentration", dimensions=dimensions, n_components=3) + concentration = Field(name="concentration", + dimensions=dimensions, n_components=3) assert_shape(concentration.shape, (11, 6, 3), "Vector field shape") - assert_equal(concentration.n_dof, 11 * 6 * 3, "Vector field degrees of freedom") + assert_equal(concentration.n_dof, 11 * 6 * 3, + "Vector field degrees of freedom") # Custom data data = np.ones((11, 6, 3)) - concentration_with_data = Field(name="concentration", dimensions=dimensions, n_components=3, data=data) - assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") + concentration_with_data = Field( + name="concentration", dimensions=dimensions, n_components=3, data=data) + assert_shape(concentration_with_data.shape, + (11, 6, 3), "Custom data field shape") + + with pytest.raises(ValueError): + viscosity.data = np.ones((1, 2, 3)) # %% Plotting @@ -41,13 +50,15 @@ def test_field_initialization(): def test_plotting(): # 1D Plot dimensions = {"x": np.linspace(0, 10, 11)} - field_1D = Field(name="1D Field", dimensions=dimensions, n_components=2, data=np.random.random((11, 2))) + field_1D = Field(name="1D Field", dimensions=dimensions, + n_components=2, data=np.random.random((11, 2))) fig, ax = field_1D.plot() assert isinstance(ax, plt.Axes), "1D plot returns one axis" # 2D Plot dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} - field_2D = Field(name="2D Field", dimensions=dimensions, n_components=3, data=np.random.random((11, 6, 3))) + field_2D = Field(name="2D Field", dimensions=dimensions, + n_components=3, data=np.random.random((11, 6, 3))) fig, axes = field_2D.plot() assert len(axes) == 3, "2D plot returns one axis per component" @@ -55,18 +66,22 @@ def test_plotting(): # %% Slicing def test_field_slicing(): - dimensions = {"axial": np.linspace(0, 10, 11), "radial": np.linspace(0, 5, 6)} + dimensions = {"axial": np.linspace( + 0, 10, 11), "radial": np.linspace(0, 5, 6)} field = Field(name="concentration", dimensions=dimensions, n_components=3) # Slice along one dimension field_sliced = field[{"axial": 0}] - assert_equal(len(field_sliced.dimensions), 1, "Field slicing reduces dimensionality") + assert_equal(len(field_sliced.dimensions), 1, + "Field slicing reduces dimensionality") assert_shape(field_sliced.shape, (6, 3), "Field slicing shape") # Slice along all dimensions field_sliced_all = field[{"axial": 0, "radial": 0}] - assert_equal(len(field_sliced_all.dimensions), 0, "Full slicing removes all dimensions") - assert_shape(field_sliced_all.shape, (3,), "Full slicing results in vector") + assert_equal(len(field_sliced_all.dimensions), 0, + "Full slicing removes all dimensions") + assert_shape(field_sliced_all.shape, (3,), + "Full slicing results in vector") # %% Normalization @@ -89,14 +104,15 @@ def test_field_normalization(): # Test 2: Verify data normalization normalized_data = normalized_field.data - assert np.isclose(np.min(normalized_data), 0.0), "Normalized data minimum is not 0." - assert np.isclose(np.max(normalized_data), 1.0), "Normalized data maximum is not 1." + assert np.isclose(np.min(normalized_data), + 0.0), "Normalized data minimum is not 0." + assert np.isclose(np.max(normalized_data), + 1.0), "Normalized data maximum is not 1." # Test 3: Ensure original field is unchanged assert np.array_equal(field.data, z), "Original field data was modified." - # %% Interpolation and Resampling def test_temperature_use_case(): @@ -136,7 +152,8 @@ def test_interpolated_field(): "radial": np.linspace(0, 5, 6) } data = np.random.random((11, 6, 3)) - concentration = Field(name="concentration", dimensions=dimensions, n_components=3, data=data) + concentration = Field(name="concentration", + dimensions=dimensions, n_components=3, data=data) # Interpolated field interp_field = FieldInterpolator(concentration) @@ -146,7 +163,8 @@ def test_interpolated_field(): def test_resampling(): dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} - field = Field(name="concentration", dimensions=dimensions, n_components=2, data=np.random.random((11, 6, 2))) + field = Field(name="concentration", dimensions=dimensions, + n_components=2, data=np.random.random((11, 6, 2))) # Resample one dimension resampled_field = field.resample({"x": 50}) @@ -154,7 +172,8 @@ def test_resampling(): # Resample all dimensions resampled_field_all = field.resample({"x": 50, "y": 25}) - assert_shape(resampled_field_all.shape, (50, 25, 2), "Resampling all dimensions") + assert_shape(resampled_field_all.shape, (50, 25, 2), + "Resampling all dimensions") def test_field_interpolation_and_derivatives(): @@ -203,6 +222,5 @@ def test_field_interpolation_and_derivatives(): # %% Run tests -import pytest if __name__ == "__main__": pytest.main('test_field.py') From ea126c710c98a3c8d0238f5182f09e3b8f14089a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Sun, 8 Dec 2024 13:27:03 +0100 Subject: [PATCH 05/13] assert n_dimensions, n_cells --- tests/test_field.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_field.py b/tests/test_field.py index 717867e..500874e 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -41,6 +41,9 @@ def test_field_initialization(): assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") + assert_equal(viscosity.n_dimensions, 3) + assert_equal(viscosity.n_cells, 11 * 6) + with pytest.raises(ValueError): viscosity.data = np.ones((1, 2, 3)) From 040c731c65f83ae1eeaa171730b113a4fdda146d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Sun, 8 Dec 2024 13:30:00 +0100 Subject: [PATCH 06/13] fix n_dimensions assertion --- tests/test_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_field.py b/tests/test_field.py index 500874e..49c5fb3 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -41,7 +41,7 @@ def test_field_initialization(): assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") - assert_equal(viscosity.n_dimensions, 3) + assert_equal(viscosity.n_dimensions, 2) assert_equal(viscosity.n_cells, 11 * 6) with pytest.raises(ValueError): From a39b16626f015ebb0b49cd21b59aafa75281d515 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Mon, 9 Dec 2024 08:46:47 +0100 Subject: [PATCH 07/13] test Field.data_flat --- tests/test_field.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/test_field.py b/tests/test_field.py index 49c5fb3..91eef5e 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -40,10 +40,17 @@ def test_field_initialization(): name="concentration", dimensions=dimensions, n_components=3, data=data) assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") + assert_equal(concentration_with_data.data_flat, np.ones(11 * 6 * 3)) assert_equal(viscosity.n_dimensions, 2) assert_equal(viscosity.n_cells, 11 * 6) + viscosity.data_flat = np.ones(11 * 6) + assert_equal(viscosity.data, np.ones((11, 6))) + + with pytest.raises(ValueError): + viscosity.data_flat = np.ones(42) + with pytest.raises(ValueError): viscosity.data = np.ones((1, 2, 3)) From 96ea7716c3616c43b8a1972af96e8beca5523e20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Mon, 9 Dec 2024 12:35:34 +0100 Subject: [PATCH 08/13] fix lint: add docstrings, remove unused code --- CADETPythonSimulator/field.py | 151 +++++++++++++++++----------------- tests/test_field.py | 115 ++++++++++++++------------ 2 files changed, 140 insertions(+), 126 deletions(-) diff --git a/CADETPythonSimulator/field.py b/CADETPythonSimulator/field.py index 8cc784b..59875df 100644 --- a/CADETPythonSimulator/field.py +++ b/CADETPythonSimulator/field.py @@ -26,12 +26,12 @@ class Field: """ def __init__( - self, - name: str, - dimensions: dict[str, npt.ArrayLike], - n_components: Optional[int] = None, - data: Optional[npt.ArrayLike] = None, - ): + self, + name: str, + dimensions: dict[str, npt.ArrayLike], + n_components: Optional[int] = None, + data: Optional[npt.ArrayLike] = None, + ): """Construct the object.""" self.name = name self.dimensions = {k: np.array(v) for k, v in dimensions.items()} @@ -107,13 +107,13 @@ def data_flat(self, data_flat: np.ndarray) -> None: self.data = data_flat.reshape(self.shape) def plot( - self, - title: str = None, - fixed_dims: Optional[dict[str, Union[float, int]]] = None, - method: str = "surface", - fig: Optional[plt.Figure] = None, - axes: Optional[plt.Axes | list[plt.Axes]] = None, - ) -> tuple[plt.Figure, plt.Axes | list[plt.Axes]]: + self, + title: str = None, + fixed_dims: Optional[dict[str, Union[float, int]]] = None, + method: str = "surface", + fig: Optional[plt.Figure] = None, + axes: Optional[plt.Axes | list[plt.Axes]] = None, + ) -> tuple[plt.Figure, plt.Axes | list[plt.Axes]]: """ Plot the field data, supporting 1D, 2D, and slicing for higher dimensions. @@ -163,7 +163,8 @@ def plot( y = self.data component_names = ( [f"Component {i}" for i in range(self.n_components)] - if self.n_components else ["Value"] + if self.n_components + else ["Value"] ) return plot_1D(x, y, component_names, title or self.name, fig, axes) @@ -173,7 +174,8 @@ def plot( z = self.data component_names = ( [f"Component {i}" for i in range(self.n_components)] - if self.n_components else ["Value"] + if self.n_components + else ["Value"] ) if method == "surface": return plot_2D_surface( @@ -186,10 +188,7 @@ def plot( else: raise ValueError(f"Unknown plotting method: {method}") - def resample( - self, - resample_dims: dict[str, Union[int, npt.ArrayLike]] - ) -> "Field": + def resample(self, resample_dims: dict[str, Union[int, npt.ArrayLike]]) -> "Field": """ Resample the field data to new dimensions using interpolation. @@ -309,10 +308,7 @@ def anti_derivative(self, dimension: str, initial_value: float = 0) -> "Field": interpolated_field = FieldInterpolator(self) return interpolated_field.anti_derivative(dimension, initial_value) - def __getitem__( - self, - slices: dict[str, float | int] - ) -> Union["Field", np.ndarray]: + def __getitem__(self, slices: dict[str, float | int]) -> Union["Field", np.ndarray]: """ Slice the field along specified dimensions. @@ -362,14 +358,18 @@ def __getitem__( n_components=self.n_components, data=sliced_data, ) + def __repr__(self) -> str: """Represantator.""" components = ( f", components={self.n_components}" - if self.n_components else ", scalar=True" + if self.n_components + else ", scalar=True" + ) + return ( + f"Field(name='{self.name}')" + + f", dimensions={list(self.dimensions.keys())}{components}" ) - return f"Field(name='{self.name}')"\ - +f", dimensions={list(self.dimensions.keys())}{components}" class NormalizedField(Field): @@ -428,7 +428,7 @@ def _initialize_interpolators(self): if self.field.n_components: self._interpolators = [ RegularGridInterpolator( - grid_points, self.field.data[..., i], method='pchip' + grid_points, self.field.data[..., i], method="pchip" ) for i in range(self.field.n_components) ] @@ -455,7 +455,7 @@ def __call__(self, **kwargs: npt.ArrayLike) -> np.ndarray: # Separate provided and missing dimensions provided_dims = set(kwargs.keys()) all_dims = set(self.field.dimensions.keys()) - missing_dims = all_dims - provided_dims + # missing_dims = all_dims - provided_dims # Validate dimensions if not provided_dims.issubset(all_dims): @@ -485,8 +485,7 @@ def __call__(self, **kwargs: npt.ArrayLike) -> np.ndarray: return result.reshape(output_shape) def _determine_output_shape( - self, - query: dict[str, Union[float, npt.ArrayLike]] + self, query: dict[str, Union[float, npt.ArrayLike]] ) -> tuple[int, ...]: """ Determine the shape of the interpolated result based on the query. @@ -543,10 +542,13 @@ def derivative(self, dimension: str): # Compute derivative for multi-component fields if self.field.n_components: - derivative_data = np.stack([ - np.gradient(self.field.data[..., i], coords, axis=axis) - for i in range(self.field.n_components) - ], axis=-1) + derivative_data = np.stack( + [ + np.gradient(self.field.data[..., i], coords, axis=axis) + for i in range(self.field.n_components) + ], + axis=-1, + ) else: derivative_data = np.gradient(self.field.data, coords, axis=axis) @@ -555,7 +557,7 @@ def derivative(self, dimension: str): name=f"{self.field.name}_derivative_{dimension}", dimensions=self.field.dimensions, n_components=self.field.n_components, - data=derivative_data + data=derivative_data, ) def anti_derivative(self, dimension: str, initial_value: float = 0): @@ -583,23 +585,26 @@ def anti_derivative(self, dimension: str, initial_value: float = 0): # Compute anti-derivative for multi-component fields if self.field.n_components: - anti_derivative_data = np.stack([ - np.cumsum( - self.field.data[..., i] * np.gradient(coords), axis=axis - ) + initial_value - for i in range(self.field.n_components) - ], axis=-1) + anti_derivative_data = np.stack( + [ + np.cumsum(self.field.data[..., i] * np.gradient(coords), axis=axis) + + initial_value + for i in range(self.field.n_components) + ], + axis=-1, + ) else: - anti_derivative_data = np.cumsum( - self.field.data * np.gradient(coords), axis=axis - ) + initial_value + anti_derivative_data = ( + np.cumsum(self.field.data * np.gradient(coords), axis=axis) + + initial_value + ) # Return new Field for the anti-derivative return Field( name=f"{self.field.name}_anti_derivative_{dimension}", dimensions=self.field.dimensions, n_components=self.field.n_components, - data=anti_derivative_data + data=anti_derivative_data, ) def __repr__(self) -> str: @@ -608,13 +613,13 @@ def __repr__(self) -> str: def plot_1D( - x: np.ndarray, - y: np.ndarray, - component_names: list[str], - title: str, - fig: Optional[plt.Figure] = None, - ax: Optional[plt.Axes] = None, - ) -> tuple[plt.Figure, plt.Axes]: + x: np.ndarray, + y: np.ndarray, + component_names: list[str], + title: str, + fig: Optional[plt.Figure] = None, + ax: Optional[plt.Axes] = None, +) -> tuple[plt.Figure, plt.Axes]: """ Plot a 1D field with a line for each component. @@ -654,14 +659,14 @@ def plot_1D( def plot_2D_surface( - x: np.ndarray, - y: np.ndarray, - z: np.ndarray, - component_names: list[str], - title: str, - fig: Optional[plt.Figure] = None, - axes: Optional[list[plt.Axes]] = None, - ) -> tuple[plt.Figure, list[plt.Axes]]: + x: np.ndarray, + y: np.ndarray, + z: np.ndarray, + component_names: list[str], + title: str, + fig: Optional[plt.Figure] = None, + axes: Optional[list[plt.Axes]] = None, +) -> tuple[plt.Figure, list[plt.Axes]]: """ Plot a 2D field with a surface for each component. @@ -692,9 +697,10 @@ def plot_2D_surface( if axes is None: fig, axes = plt.subplots( - 1, n_components, + 1, + n_components, figsize=(5 * n_components, 5), - subplot_kw={'projection': '3d'} + subplot_kw={"projection": "3d"}, ) axes = [axes] if n_components == 1 else axes @@ -710,14 +716,14 @@ def plot_2D_surface( def plot_2D_heatmap( - x: np.ndarray, - y: np.ndarray, - z: np.ndarray, - component_names: list[str], - title: str, - fig: Optional[plt.Figure] = None, - axes: Optional[list[plt.Axes]] = None, - ) -> tuple[plt.Figure, list[plt.Axes]]: + x: np.ndarray, + y: np.ndarray, + z: np.ndarray, + component_names: list[str], + title: str, + fig: Optional[plt.Figure] = None, + axes: Optional[list[plt.Axes]] = None, +) -> tuple[plt.Figure, list[plt.Axes]]: """ Plot a 2D field with a heatmap for each component. @@ -748,9 +754,7 @@ def plot_2D_heatmap( if axes is None: fig, axes = plt.subplots( - 1, n_components, - figsize=(5 * n_components, 5), - constrained_layout=True + 1, n_components, figsize=(5 * n_components, 5), constrained_layout=True ) axes = [axes] if n_components == 1 else axes @@ -763,4 +767,3 @@ def plot_2D_heatmap( ax.set_ylabel("Y-axis") return fig, axes - diff --git a/tests/test_field.py b/tests/test_field.py index 91eef5e..adedec3 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -6,40 +6,39 @@ def assert_equal(value, expected, message=""): + """Assert equality.""" message = f"Test failed: {message}. Expected {expected}, got {value}." assert value == expected, message def assert_shape(shape, expected_shape, message=""): + """Assert equality of shapes.""" message = f"Test failed: {message}. Expected shape {expected_shape}, got {shape}." assert shape == expected_shape, message # %% Initialization + def test_field_initialization(): + """Test initialization of Field class.""" # Scalar field - dimensions = { - "axial": np.linspace(0, 10, 11), - "radial": np.linspace(0, 5, 6) - } + dimensions = {"axial": np.linspace(0, 10, 11), "radial": np.linspace(0, 5, 6)} viscosity = Field(name="viscosity", dimensions=dimensions) assert_shape(viscosity.shape, (11, 6), "Scalar field shape") assert_equal(viscosity.n_dof, 11 * 6, "Scalar field degrees of freedom") # Vector field - concentration = Field(name="concentration", - dimensions=dimensions, n_components=3) + concentration = Field(name="concentration", dimensions=dimensions, n_components=3) assert_shape(concentration.shape, (11, 6, 3), "Vector field shape") - assert_equal(concentration.n_dof, 11 * 6 * 3, - "Vector field degrees of freedom") + assert_equal(concentration.n_dof, 11 * 6 * 3, "Vector field degrees of freedom") # Custom data data = np.ones((11, 6, 3)) concentration_with_data = Field( - name="concentration", dimensions=dimensions, n_components=3, data=data) - assert_shape(concentration_with_data.shape, - (11, 6, 3), "Custom data field shape") + name="concentration", dimensions=dimensions, n_components=3, data=data + ) + assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") assert_equal(concentration_with_data.data_flat, np.ones(11 * 6 * 3)) assert_equal(viscosity.n_dimensions, 2) @@ -57,46 +56,60 @@ def test_field_initialization(): # %% Plotting + def test_plotting(): + """Test Field plotting results.""" # 1D Plot dimensions = {"x": np.linspace(0, 10, 11)} - field_1D = Field(name="1D Field", dimensions=dimensions, - n_components=2, data=np.random.random((11, 2))) + field_1D = Field( + name="1D Field", + dimensions=dimensions, + n_components=2, + data=np.random.random((11, 2)), + ) fig, ax = field_1D.plot() assert isinstance(ax, plt.Axes), "1D plot returns one axis" # 2D Plot dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} - field_2D = Field(name="2D Field", dimensions=dimensions, - n_components=3, data=np.random.random((11, 6, 3))) + field_2D = Field( + name="2D Field", + dimensions=dimensions, + n_components=3, + data=np.random.random((11, 6, 3)), + ) fig, axes = field_2D.plot() assert len(axes) == 3, "2D plot returns one axis per component" # %% Slicing + def test_field_slicing(): - dimensions = {"axial": np.linspace( - 0, 10, 11), "radial": np.linspace(0, 5, 6)} + """Test Field slicing.""" + dimensions = {"axial": np.linspace(0, 10, 11), "radial": np.linspace(0, 5, 6)} field = Field(name="concentration", dimensions=dimensions, n_components=3) # Slice along one dimension field_sliced = field[{"axial": 0}] - assert_equal(len(field_sliced.dimensions), 1, - "Field slicing reduces dimensionality") + assert_equal( + len(field_sliced.dimensions), 1, "Field slicing reduces dimensionality" + ) assert_shape(field_sliced.shape, (6, 3), "Field slicing shape") # Slice along all dimensions field_sliced_all = field[{"axial": 0, "radial": 0}] - assert_equal(len(field_sliced_all.dimensions), 0, - "Full slicing removes all dimensions") - assert_shape(field_sliced_all.shape, (3,), - "Full slicing results in vector") + assert_equal( + len(field_sliced_all.dimensions), 0, "Full slicing removes all dimensions" + ) + assert_shape(field_sliced_all.shape, (3,), "Full slicing results in vector") # %% Normalization + def test_field_normalization(): + """Test Field normalization.""" # Define dimensions and data x = np.linspace(0, 10, 11) y = np.linspace(0, 5, 6) @@ -114,10 +127,8 @@ def test_field_normalization(): # Test 2: Verify data normalization normalized_data = normalized_field.data - assert np.isclose(np.min(normalized_data), - 0.0), "Normalized data minimum is not 0." - assert np.isclose(np.max(normalized_data), - 1.0), "Normalized data maximum is not 1." + assert np.isclose(np.min(normalized_data), 0.0), "Normalized data minimum is not 0." + assert np.isclose(np.max(normalized_data), 1.0), "Normalized data maximum is not 1." # Test 3: Ensure original field is unchanged assert np.array_equal(field.data, z), "Original field data was modified." @@ -125,7 +136,9 @@ def test_field_normalization(): # %% Interpolation and Resampling + def test_temperature_use_case(): + """Test examplary use case of interpolating temperature from a Field.""" temperature_profile_dimensions = { "time": np.linspace(0, 3600, 3601), "axial": np.linspace(0, 0.5, 6), @@ -148,22 +161,14 @@ def calculate_temperature_at_t_x(t, x): def calculate_adsorption_from_temperature(k_0, k_1, T): return k_0 * np.exp(k_1 / T) - def calculate_parameter(t, x): - T = calculate_temperature_at_t_x(t, x) - k = calculate_adsorption_from_temperature(k_0, k_1, T) - - def model_equation(t): - ... - def test_interpolated_field(): - dimensions = { - "axial": np.linspace(0, 10, 11), - "radial": np.linspace(0, 5, 6) - } + """Test interpolation of a field.""" + dimensions = {"axial": np.linspace(0, 10, 11), "radial": np.linspace(0, 5, 6)} data = np.random.random((11, 6, 3)) - concentration = Field(name="concentration", - dimensions=dimensions, n_components=3, data=data) + concentration = Field( + name="concentration", dimensions=dimensions, n_components=3, data=data + ) # Interpolated field interp_field = FieldInterpolator(concentration) @@ -172,9 +177,14 @@ def test_interpolated_field(): def test_resampling(): + """Test field resampling.""" dimensions = {"x": np.linspace(0, 10, 11), "y": np.linspace(0, 5, 6)} - field = Field(name="concentration", dimensions=dimensions, - n_components=2, data=np.random.random((11, 6, 2))) + field = Field( + name="concentration", + dimensions=dimensions, + n_components=2, + data=np.random.random((11, 6, 2)), + ) # Resample one dimension resampled_field = field.resample({"x": 50}) @@ -182,14 +192,14 @@ def test_resampling(): # Resample all dimensions resampled_field_all = field.resample({"x": 50, "y": 25}) - assert_shape(resampled_field_all.shape, (50, 25, 2), - "Resampling all dimensions") + assert_shape(resampled_field_all.shape, (50, 25, 2), "Resampling all dimensions") def test_field_interpolation_and_derivatives(): + """Test field interpolation and derivatives.""" # Define test dimensions and data x = np.linspace(0, 10, 11) # 11 points along x - y = np.linspace(0, 5, 6) # 6 points along y + y = np.linspace(0, 5, 6) # 6 points along y z = np.outer(np.sin(x), np.cos(y)) # Generate 2D scalar field data # Create a Field @@ -201,13 +211,13 @@ def test_field_interpolation_and_derivatives(): # Test 1: Evaluate at an arbitrary point eval_result = interp_field(x=2.5) assert_shape( - eval_result.shape, (6, ), - "Evaluation result should be a scalar for scalar fields." + eval_result.shape, + (6,), + "Evaluation result should be a scalar for scalar fields.", ) eval_result = interp_field(x=2.5, y=1.1) assert_shape( - eval_result.shape, (), - "Evaluation result should be a scalar for scalar fields." + eval_result.shape, (), "Evaluation result should be a scalar for scalar fields." ) # Test 2: Compute derivative along x @@ -225,12 +235,13 @@ def test_field_interpolation_and_derivatives(): # Test 5: Edge Case - Evaluate at boundary boundary_value = interp_field(x=0, y=5) assert_shape( - boundary_value.shape, (), - "Boundary evaluation should return a scalar for scalar fields." + boundary_value.shape, + (), + "Boundary evaluation should return a scalar for scalar fields.", ) # %% Run tests if __name__ == "__main__": - pytest.main('test_field.py') + pytest.main("test_field.py") From 8585a8ff85c7a7c8a947825ac51e3e9dde6cb02e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Mon, 9 Dec 2024 15:23:16 +0100 Subject: [PATCH 09/13] use assert_array_equal for arrays --- tests/test_field.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/test_field.py b/tests/test_field.py index adedec3..3a5c86f 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -2,6 +2,7 @@ import numpy as np from CADETPythonSimulator.field import Field, FieldInterpolator import matplotlib.pyplot as plt +from numpy.testing import assert_array_equal # %% Testing utilities @@ -39,13 +40,13 @@ def test_field_initialization(): name="concentration", dimensions=dimensions, n_components=3, data=data ) assert_shape(concentration_with_data.shape, (11, 6, 3), "Custom data field shape") - assert_equal(concentration_with_data.data_flat, np.ones(11 * 6 * 3)) + assert_array_equal(concentration_with_data.data_flat, np.ones(11 * 6 * 3)) assert_equal(viscosity.n_dimensions, 2) assert_equal(viscosity.n_cells, 11 * 6) viscosity.data_flat = np.ones(11 * 6) - assert_equal(viscosity.data, np.ones((11, 6))) + assert_array_equal(viscosity.data, np.ones((11, 6))) with pytest.raises(ValueError): viscosity.data_flat = np.ones(42) From 03f4cbe15c399fe98f27d5689d82a22dd1895b40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Sat, 14 Dec 2024 23:55:36 +0100 Subject: [PATCH 10/13] cover 3d plot --- tests/test_field.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/test_field.py b/tests/test_field.py index 3a5c86f..9e29af0 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -82,6 +82,24 @@ def test_plotting(): fig, axes = field_2D.plot() assert len(axes) == 3, "2D plot returns one axis per component" + # plot 3D field + dimensions = { + "x": np.linspace(0, 10, 11), + "y": np.linspace(0, 5, 6), + "z": np.linspace(0, 2, 3), + } + field_3D = Field( + name="3D Field", + dimensions=dimensions, + n_components=4, + data=np.random.random((11, 6, 3, 4)), + ) + with pytest.raises(ValueError): + field_3D.plot() + + fig, axis = field_3D.plot(fixed_dims={"x": 1, "y": 1}) + assert isinstance(axis, plt.Axes), "3D plot with two fixed dimensions has one axis" + # %% Slicing From 53ae102376a190f3e1f6571e37e620d3820009f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Sun, 15 Dec 2024 16:37:03 +0100 Subject: [PATCH 11/13] read csv --- tests/test_field.py | 13 +++++++++++++ tests/test_temp.csv | 11 +++++++++++ 2 files changed, 24 insertions(+) create mode 100644 tests/test_temp.csv diff --git a/tests/test_field.py b/tests/test_field.py index 9e29af0..a78575e 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -260,6 +260,19 @@ def test_field_interpolation_and_derivatives(): ) +def test_field_from_csv(): + """Test reading field from csv.""" + data = np.genfromtxt("test_temp.csv", delimiter=",", names=True) + dim_time = np.unique(data["time"]) + dim_axial = np.unique(data["axial"]) + values = data["Temperature"].reshape(len(dim_time), len(dim_axial)) + field = Field( + name="temp_field", + dimensions={"time": dim_time, "axial": dim_axial}, + data=values, + ) + + # %% Run tests if __name__ == "__main__": diff --git a/tests/test_temp.csv b/tests/test_temp.csv new file mode 100644 index 0000000..f672bd9 --- /dev/null +++ b/tests/test_temp.csv @@ -0,0 +1,11 @@ +"time", "axial", "Temperature" +0, 0.0, 20 +0, 0.5, 20 +0, 1.0, 20 +1800, 0.0, 25 +1800, 0.5, 25 +1800, 1.0, 25 +3600, 0.0, 30 +3600, 0.5, 30 +3600, 1.0, 30 + From 0b93c5db47e4a7a44c0caa5589b4df5c6266ad61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Mon, 16 Dec 2024 10:17:46 +0100 Subject: [PATCH 12/13] fix path --- tests/test_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_field.py b/tests/test_field.py index a78575e..72abb9d 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -262,7 +262,7 @@ def test_field_interpolation_and_derivatives(): def test_field_from_csv(): """Test reading field from csv.""" - data = np.genfromtxt("test_temp.csv", delimiter=",", names=True) + data = np.genfromtxt("tests/test_temp.csv", delimiter=",", names=True) dim_time = np.unique(data["time"]) dim_axial = np.unique(data["axial"]) values = data["Temperature"].reshape(len(dim_time), len(dim_axial)) From 366a12780dceb27a265590c2521d978b1c788649 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20H=C3=BClsmann?= Date: Mon, 16 Dec 2024 17:25:50 +0100 Subject: [PATCH 13/13] shuffle csv, values by index --- tests/test_field.py | 20 +++++++++++++++++--- tests/test_temp.csv | 24 +++++++++++++----------- 2 files changed, 30 insertions(+), 14 deletions(-) diff --git a/tests/test_field.py b/tests/test_field.py index 72abb9d..945c7e2 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -263,14 +263,28 @@ def test_field_interpolation_and_derivatives(): def test_field_from_csv(): """Test reading field from csv.""" data = np.genfromtxt("tests/test_temp.csv", delimiter=",", names=True) - dim_time = np.unique(data["time"]) - dim_axial = np.unique(data["axial"]) - values = data["Temperature"].reshape(len(dim_time), len(dim_axial)) + + # Initialize dimensions + dim_time, idx_time = np.unique(data["time"], return_inverse=True) + dim_axial, idx_axial = np.unique(data["axial"], return_inverse=True) + + # Initialize values + values = np.zeros((len(dim_time), len(dim_axial))) + for idx in range(len(data)): + values[idx_time[idx]][idx_axial[idx]] = data[idx]["Temperature"] + field = Field( name="temp_field", dimensions={"time": dim_time, "axial": dim_axial}, data=values, ) + assert_shape(field.shape, (4, 3), "Field shape mismatch") + + interp_field = FieldInterpolator(field) + assert_equal(interp_field(time=1700, axial=0), 24, "Wrong interpolated temperature") + assert_equal( + interp_field(time=300, axial=0.3), 17.34, "Wrong interpolated temperature" + ) # %% Run tests diff --git a/tests/test_temp.csv b/tests/test_temp.csv index f672bd9..a9ac3d3 100644 --- a/tests/test_temp.csv +++ b/tests/test_temp.csv @@ -1,11 +1,13 @@ -"time", "axial", "Temperature" -0, 0.0, 20 -0, 0.5, 20 -0, 1.0, 20 -1800, 0.0, 25 -1800, 0.5, 25 -1800, 1.0, 25 -3600, 0.0, 30 -3600, 0.5, 30 -3600, 1.0, 30 - +"time","axial","Temperature" +0,0.0,20 +0,0.5,15 +1500,0.0,22 +1500,0.5,16.5 +1800,0.0,25 +1500,1.0,11 +3600,0.5,22.5 +0,1.0,10 +1800,1.0,12.5 +3600,0.0,30 +1800,0.5,18.75 +3600,1.0,15