Source code for openmc.data.grid

import numpy as np


[docs]def linearize(x, f, tolerance=0.001): """Return a tabulated representation of a one-variable function Parameters ---------- x : Iterable of float Initial x values at which the function should be evaluated f : Callable Function of a single variable tolerance : float Tolerance on the interpolation error Returns ------- numpy.ndarray Tabulated values of the independent variable numpy.ndarray Tabulated values of the dependent variable """ # Make sure x is a numpy array x = np.asarray(x) # Initialize output arrays x_out = [] y_out = [] # Initialize stack x_stack = [x[0]] y_stack = [f(x[0])] for i in range(x.shape[0] - 1): x_stack.insert(0, x[i + 1]) y_stack.insert(0, f(x[i + 1])) while True: # Get the bounding points currently on the stack x_high, x_low = x_stack[-2:] y_high, y_low = y_stack[-2:] # Evaluate the function at the midpoint x_mid = 0.5*(x_low + x_high) y_mid = f(x_mid) # Linearly interpolate between the bounding points y_interp = y_low + (y_high - y_low)/(x_high - x_low)*(x_mid - x_low) # Check the error on the interpolated point and compare to tolerance error = abs((y_interp - y_mid)/y_mid) if error > tolerance: x_stack.insert(-1, x_mid) y_stack.insert(-1, y_mid) else: x_out.append(x_stack.pop()) y_out.append(y_stack.pop()) if len(x_stack) == 1: break x_out.append(x_stack.pop()) y_out.append(y_stack.pop()) return np.array(x_out), np.array(y_out)
[docs]def thin(x, y, tolerance=0.001): """Check for (x,y) points that can be removed. Parameters ---------- x : numpy.ndarray Independent variable y : numpy.ndarray Dependent variable tolerance : float Tolerance on interpolation error Returns ------- numpy.ndarray Tabulated values of the independent variable numpy.ndarray Tabulated values of the dependent variable """ # Initialize output arrays x_out = x.copy() y_out = y.copy() N = x.shape[0] i_left = 0 i_right = 2 while i_left < N - 2 and i_right < N: m = (y[i_right] - y[i_left])/(x[i_right] - x[i_left]) for i in range(i_left + 1, i_right): # Determine error in interpolated point y_interp = y[i_left] + m*(x[i] - x[i_left]) if abs(y[i]) > 0.: error = abs((y_interp - y[i])/y[i]) else: error = 2*tolerance if error > tolerance: for i_remove in range(i_left + 1, i_right - 1): x_out[i_remove] = np.nan y_out[i_remove] = np.nan i_left = i_right - 1 i_right = i_left + 1 break i_right += 1 for i_remove in range(i_left + 1, i_right - 1): x_out[i_remove] = np.nan y_out[i_remove] = np.nan return x_out[np.isfinite(x_out)], y_out[np.isfinite(y_out)]