New tabulation framework

Common/CommonMethods.py

@@ -203,4 +203,73 @@ def writeMLPForSU2(file_out:str, weights:list[np.ndarray], biases:list[np.ndarra
         fid.write("\t".join("%+.16e" % float(b) for b in B) + "\n")
 
     fid.close()
-    return
+    return
+
+def shoelace(XY:np.ndarray[float]):
+    """Shoelace algorithm for area computations
+
+    :param XY: hull node coordinates
+    :type XY: np.ndarray[float]
+    :return: area of concave hull
+    :rtype: float
+    """
+    x = XY[:,0]
+    y = XY[:,1]
+    S1 = np.sum(x*np.roll(y,-1))
+    S2 = np.sum(y*np.roll(x,-1))
+
+    area = .5*np.absolute(S1 - S2)
+    return area
+
+
+def FiniteDifferenceDerivative(y:np.ndarray[float], x:np.ndarray[float]):
+    """Calculate second-order accurate, one-dimensional finite-difference derivatives of y with respect to x.
+
+    :param y: data to calculate the finite-differences for.
+    :type y: np.ndarray[float]
+    :param x: axial coordinates.
+    :type x: np.ndarray[float]
+    :return: finite-difference derivatives of y with respect to x.
+    :rtype: np.ndarray[float]
+    """
+    Np = len(x)
+    dydx = np.zeros(Np)
+    for i in range(1, Np-1):
+        y_m = y[i-1]
+        y_p = y[i+1]
+        y_0 = y[i]
+        x_m = x[i-1]
+        x_p = x[i+1]
+        x_0 = x[i]
+        dx_1 = x_p - x_0
+        dx_2 = x_0 - x_m
+        dx2_1 = dx_1*dx_1
+        dx2_2 = dx_2*dx_2
+        if (dx_1==0) or (dx_2==0):
+            dydx[i] = 0.0
+        else:
+            dydx[i] = (dx2_2 * y_p + (dx2_1 - dx2_2)*y_0 - dx2_1*y_m)/(dx_1*dx_2*(dx_1+dx_2))
+    dx_1 = x[1] - x[0]
+    dx_2 = x[2] - x[0]
+    dx2_1 = dx_1*dx_1
+    dx2_2 = dx_2*dx_2
+    y_0 = y[0]
+    y_p = y[1]
+    y_pp = y[2]
+    if (dx_1==0) or (dx_2==0):
+        dydx[0] = 0.0
+    else:
+        dydx[0] = (dx2_1 * y_pp + (dx2_2 - dx2_1)*y_0 - dx2_2*y_p)/(dx_1*dx_2*(dx_1 - dx_2))
+
+    dx_1 = x[-2] - x[-1]
+    dx_2 = x[-3] - x[-1]
+    dx2_1 = dx_1*dx_1
+    dx2_2 = dx_2*dx_2
+    y_0 = y[-1]
+    y_p = y[-2]
+    y_pp = y[-3]
+    if (dx_1==0) or (dx_2==0):
+        dydx[-1] = 0.0
+    else:
+        dydx[-1] = (dx2_1 * y_pp + (dx2_2 - dx2_1)*y_0 - dx2_2*y_p)/(dx_1*dx_2*(dx_1 - dx_2))
+    return dydx

Common/DataDrivenConfig.py

@@ -38,7 +38,7 @@
 #---------------------------------------------------------------------------------------------#
 # Importing DataMiner classes and functions
 #---------------------------------------------------------------------------------------------#
-from Common.Properties import DefaultSettings_NICFD, DefaultSettings_FGM, FlameletSolverOptions
+from Common.Properties import DefaultSettings_NICFD, DefaultSettings_FGM, FGMVars, FlameletSolverOptions
 from Common.Config_base import Config
 from Common.CommonMethods import *
 
@@ -925,6 +925,10 @@ def __SynchronizeSettings(self):
 
         self.SetAverageLewisNumbers()
 
+        if self.__mix_status_lower==self.__mix_status_upper:
+            self.SetControllingVariables([FGMVars.ProgressVariable.name, FGMVars.EnthalpyTot.name])
+        else:
+            self.SetControllingVariables([FGMVars.ProgressVariable.name, FGMVars.EnthalpyTot.name, FGMVars.MixtureFraction.name])
         return
 
     def GetConstSpecieLewisNumbers(self):

Common/Interpolators.py

@@ -25,7 +25,9 @@
 
 from scipy.spatial import cKDTree as KDTree
 import numpy as np
-
+import pandas as pd
+from tqdm import tqdm
+from sklearn.preprocessing import MinMaxScaler
 class Invdisttree:
     """ inverse-distance-weighted interpolation using KDTree:
 invdisttree = Invdisttree( X, z )  -- data points, values
@@ -77,41 +79,251 @@ class Invdisttree:
     """
 # anykernel( dj / av dj ) is also scale-free
 # error analysis, |f(x) - idw(x)| ? todo: regular grid, nnear ndim+1, 2*ndim
+    __state_space_data:np.ndarray[float]=None
+    __nDim_input:int=None
+    __number_of_states:int=None
+    __KD_tree:KDTree = None
+    __single_query_sample:bool=False
+    __number_of_nearest_neighbors:int=None
+    __inverse_distance_exponent:float=None
+
+    def __init__( self, samples_state_space:np.ndarray[float], state_data:np.ndarray[float], leafsize:int=10):
+        assert len(samples_state_space) == len(state_data), "len(X) %d != len(z) %d" % (len(samples_state_space), len(state_data))
+        self.__nDim_input = np.shape(samples_state_space)[1]
+        self.__number_of_states = np.shape(state_data)[1]
+        self.__KD_tree = KDTree( samples_state_space, leafsize=leafsize )
+        self.__state_space_data = state_data
+
+
+    def query( self, query_samples:np.ndarray[float], nearest_neighbors:int=3, inverse_distance_exponent:float=2 ):
+        """Interpolate based based on the inverse-distance-weighted values of a specified number of nearest neighbors.
+
+        :param query_samples: query samples
+        :type query_samples: np.ndarray[float]
+        :param nearest_neighbors: number of nearest neighbors, defaults to 3
+        :type nearest_neighbors: int, optional
+        :param inverse_distance_exponent: exponent of the inverse distance to the nearest neighbors, defaults to 2
+        :type inverse_distance_exponent: float, optional
+        :return: interpolated data from nearest neighbors
+        :rtype: np.ndarray[float]
+        """
+        self.__number_of_nearest_neighbors = nearest_neighbors
+        self.__inverse_distance_exponent = inverse_distance_exponent
+
+        formatted_query_samples = self.__parse_query_input(query_samples)
+
+        return self.__interpolate_from_inverse_distance(formatted_query_samples)
+
+    def __parse_query_input(self, input_query_samples:np.ndarray[float]):
+        formatted_query_samples = self.__format_query_input(input_query_samples)
+        self.__check_correct_input_dimensions(formatted_query_samples)
+        return formatted_query_samples
+
+    def __format_query_input(self, unformatted_query_samples:np.ndarray[float]):
+        number_of_dimensions = unformatted_query_samples.ndim
+        if number_of_dimensions == 1:
+            self.__single_query_sample = True
+            return unformatted_query_samples.reshape(1, -1)
+        else:
+            self.__single_query_sample = False
+            return unformatted_query_samples
+
+    def __check_correct_input_dimensions(self, query_samples:np.ndarray[float]):
+        if query_samples.shape[1] != self.__nDim_input:
+            raise Exception("Incorrect input dimensionality")
+        return
+
+    def __interpolate_from_inverse_distance(self, query_input_samples:np.ndarray[float]):
+        distances, nearest_neighbor_indices = self.__evaluate_KD_tree(query_input_samples, self.__number_of_nearest_neighbors)
+        interpolated_data = self.__interpolate_from_nearest_neighbors(nearest_neighbor_indices, distances)
+        if self.__single_query_sample:
+            return interpolated_data[0]
+        else:
+            return interpolated_data
+
+    def queryJacobian(self, query_input_samples:np.ndarray[float], nearest_neighbors:int=3, inverse_distance_exponent:float=2):
+
+        dx = 1e-2
+        grads = np.zeros([self.__nDim_input, query_input_samples.shape[0], self.__state_space_data.shape[1]])
+        for i in range(self.__nDim_input):
+            q_copy = query_input_samples.copy()
+            if query_input_samples.ndim==1:
+                q_copy[i] += dx
+            else:
+                q_copy[:, i] += dx
+            interp_data_plus = self.query(q_copy, nearest_neighbors, inverse_distance_exponent)
+            if query_input_samples.ndim==1:
+                q_copy[i] -= 2*dx
+            else:
+                q_copy[:, i] -= 2*dx
+            interp_data_minus = self.query(q_copy, nearest_neighbors, inverse_distance_exponent)
+
+            grads[i] = (interp_data_plus - interp_data_minus) / (2*dx)
+        return grads
+
+    def __evaluate_KD_tree(self, query_input_samples:np.ndarray[float], number_of_nearest_neighbors:int):
+        distances_to_nearest_neighbors, nearest_neighbor_indices = self.__KD_tree.query(query_input_samples, k=number_of_nearest_neighbors,eps=0)
+        return distances_to_nearest_neighbors, nearest_neighbor_indices
+
+    def __interpolate_from_nearest_neighbors(self, nearest_neighbor_indices:np.ndarray[int], distances_to_neighbors:np.ndarray[float]):
+        number_of_query_samples = np.shape(distances_to_neighbors)[0]
+        interpolated_state_data = np.zeros([number_of_query_samples, self.__number_of_states])
+        for j, (dist, ix) in enumerate(zip( distances_to_neighbors, nearest_neighbor_indices )):
+            if self.__single_nearest_neighbor():
+                interpolated_state_data[j] = self.__state_space_data[ix]
+            elif self.__nearest_sample_is_too_close(dist):
+                interpolated_state_data[j] = self.__state_space_data[ix[0]]
+            else:  # weight z s by 1/dist --
+                nearest_neighbor_coefficients = np.power(dist, -self.__inverse_distance_exponent)
+                nearest_neighbor_coefficients /= np.sum(nearest_neighbor_coefficients)
+                interpolated_state_data[j] = np.dot( nearest_neighbor_coefficients, self.__state_space_data[ix])
 
-    def __init__( self, X, z, leafsize=10, stat=0 ):
-        assert len(X) == len(z), "len(X) %d != len(z) %d" % (len(X), len(z))
-        self.tree = KDTree( X, leafsize=leafsize )  # build the tree
-        self.z = z
-        self.stat = stat
-        self.wn = 0
-        self.wsum = None
+        return interpolated_state_data
+
+    def __single_nearest_neighbor(self):
+        return self.__number_of_nearest_neighbors==1
+
+    def __nearest_sample_is_too_close(self, ordered_distances:np.ndarray[float]):
+        return ordered_distances[0] < 1e-10
+
+    def query_with_local_parameters( self, q:np.ndarray[float], nnear_vals:np.ndarray[int], p_vals:np.ndarray[float]):
+        """Interpolate with number of nearest neighbors and p-factor for each sample.
 
-    def __call__( self, q, nnear=3, eps=0, p=2, weights=None ):
+        :param q: query samples
+        :type q: np.ndarray[float]
+        :param nnear_vals: number of nearest neighbors for each query sample
+        :type nnear_vals: np.array[int]
+        :param p_vals: distance exponent factor for each query sample
+        :type p_vals: np.array[float]
+        :return: interpolated data from nearest neighbors
+        :rtype: np.ndarray[float]
+        """
             # nnear nearest neighbours of each query point --
         q = np.asarray(q)
+        max_nnear = max(max(nnear_vals), 2)
         qdim = q.ndim
         if qdim == 1:
-            q = np.array([q])
-        if self.wsum is None:
-            self.wsum = np.zeros(nnear)
+            q = np.expand_dims(q, -1)
+            nnear_vals = np.expand_dims(nnear_vals, -1)
+            p_vals = np.expand_dims(p_vals, -1)
 
-        self.distances, self.ix = self.tree.query( q, k=nnear, eps=eps )
-        interpol = np.zeros( (len(self.distances),) + np.shape(self.z[0]) )
+        eps=0
+        self.distances, self.ix = self.__KD_tree.query( q, k=max_nnear, eps=eps )
+        interpol = np.zeros([np.shape(q)[0], np.shape(self.z)[1]])
         jinterpol = 0
-        for dist, ix in zip( self.distances, self.ix ):
-            if nnear == 1:
-                wz = self.z[ix]
-            elif dist[0] < 1e-10:
+        for p, nnear, dist, ix in zip(p_vals, nnear_vals, self.distances, self.ix ):
+            if nnear == 1 or dist[0] < 1e-10:
                 wz = self.z[ix[0]]
             else:  # weight z s by 1/dist --
-                w = 1 / dist**p
-                if weights is not None:
-                    w *= weights[ix]  # >= 0
+                w = np.power(dist[:nnear], -p)
                 w /= np.sum(w)
-                wz = np.dot( w, self.z[ix])
-                if self.stat:
-                    self.wn += 1
-                    self.wsum += w
-            interpol[jinterpol] = wz
+                wz = np.dot(w,self.z[ix[:nnear]])
+            interpol[jinterpol, :] = wz
             jinterpol += 1
-        return interpol if qdim > 1  else interpol[0]
+        return interpol if qdim > 1  else interpol[0]
+
+class fluidDataInterpolator:
+
+    __controlling_variable_nodes:np.ndarray[float] = None
+    __lookup_data:np.ndarray[float] = None
+    __lookup_tree:Invdisttree = None
+    __n_near_single:int = 6
+    __p_fac_single:int = 2
+    __state_vars:list[str] = None
+
+    def __init__(self,cv_data:np.ndarray[float], state_data:pd.DataFrame, number_of_nearest_neighbors:int=None, inverse_distance_exponent:float=None):
+        """Inverse distance weighted interpolation algorithm for the interpolation of thermodynamic state data from unstructured point clouds.
+
+        :param cv_data: scaled controlling variable data from point cloud samples
+        :type cv_data: np.ndarray[float]
+        :param state_data: interpolation data from point cloud samples.
+        :type state_data: pd.DataFrame
+        :param number_of_nearest_neighbors: number of nearest neighbors from which fluid data is interpolated, defaults to None
+        :type number_of_nearest_neighbors: int, optional
+        :param inverse_distance_exponent: exponent applied to the inverse distance of the nearest neighbors, defaults to None
+        :type inverse_distance_exponent: float, optional
+        """
+        self.__controlling_variable_nodes = cv_data
+        self.__lookup_data = state_data.values
+        self.__state_vars = list(state_data.keys())
+        self.__lookup_tree = Invdisttree(cv_data, state_data.values)
+        self.__n_near_single = number_of_nearest_neighbors
+        self.__p_fac_single = inverse_distance_exponent
+
+        if not self.__n_near_single or not self.__p_fac_single:
+            self.tuneTreeParameters()
+        return
+
+    def tuneTreeParameters(self):
+
+        Np_total = np.shape(self.__controlling_variable_nodes)[0]
+        Np_train = int(0.8*Np_total)
+
+        n_control_vars = np.shape(self.__controlling_variable_nodes)[1]
+        scaler = MinMaxScaler()
+        lookup_data_scaled = scaler.fit_transform(self.__lookup_data)
+        full_data_set = np.column_stack((self.__controlling_variable_nodes, lookup_data_scaled))
+        np.random.shuffle(full_data_set)
+
+        train_data = full_data_set[:Np_train]
+        test_data = full_data_set[Np_train:]
+
+        cv_data_train= train_data[:, :n_control_vars]
+        cv_data_test= test_data[:, :n_control_vars]
+
+        train_tree = Invdisttree(cv_data_train, cv_data_train)
+        n_near_range = range(1, 20)
+        p_range = np.linspace(0, 6, 10)
+        RMS_ppv = np.zeros([len(n_near_range), len(p_range)])
+        for i in tqdm(range(len(n_near_range)),desc="Searching for optimal tree parameters..."):
+            for j in range(len(p_range)):
+                interpolated_state = train_tree.query(cv_data_test, nearest_neighbors=n_near_range[i], inverse_distance_exponent=p_range[j])
+                rms_local = np.sum(np.power(interpolated_state - cv_data_test, 2))
+                RMS_ppv[i,j] = rms_local
+        [imin,jmin] = divmod(RMS_ppv.argmin(), RMS_ppv.shape[1])
+        self.__n_near_single = n_near_range[imin]
+        self.__p_fac_single = p_range[jmin]
+
+        print("Number of nearest neighbors: %i" % self.__n_near_single)
+        print("Inverse distance exponent: %.2f" % self.__p_fac_single)
+        return
+
+    def __call__(self, cv_data_query:np.ndarray[float],varnames:list[str]=[]):
+        ndim_input = cv_data_query.ndim
+        if ndim_input < 2:
+            cv_data_query = cv_data_query[np.newaxis]
+
+        fluid_data_interp = self.__lookup_tree.query(cv_data_query, nearest_neighbors=self.__n_near_single, inverse_distance_exponent=self.__p_fac_single)
+        nvars =len(varnames)
+        if nvars > 0:
+            var_indices = [self.__state_vars.index(var) for var in varnames]
+
+            if fluid_data_interp.shape[0] == 1:
+                return fluid_data_interp[0,var_indices]
+            else:
+                return fluid_data_interp[:, var_indices]
+        else:
+            if fluid_data_interp.shape[0] == 1:
+                return fluid_data_interp[0]
+            else:
+                return fluid_data_interp
+
+    def Jacobian(self, cv_data_query:np.ndarray[float],varnames:list[str]=[]):
+        ndim_input = cv_data_query.ndim
+        if ndim_input < 2:
+            cv_data_query = cv_data_query[np.newaxis]
+
+        grads = self.__lookup_tree.queryJacobian(cv_data_query, nearest_neighbors=self.__n_near_single, inverse_distance_exponent=self.__p_fac_single)
+        nvars =len(varnames)
+        if nvars > 0:
+            var_indices = [self.__state_vars.index(var) for var in varnames]
+
+            if grads.shape[1] == 1:
+                return grads[:,0,var_indices]
+            else:
+                return grads[:,:,var_indices]
+        else:
+            if grads.shape[1] == 1:
+                return grads[:,0,:]
+            else:
+                return grads