<jupyter_start><jupyter_text>DescriptionThis notebook demonstrates how to use the dataset. A large set of features and labels is formed. A tensor basis neural network (TBNN) is created using Keras. The network is trained and evaluated on some sample points.  Loading the dataset features and labels, creating a training/validation dataset<jupyter_code>import numpy as np
np.load(f"/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1100_I1.npy", allow_pickle=True) 
np.load('/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1100_q.npy', allow_pickle=True)
import numpy as np

x = np.arange(10).reshape(1, 5, 2)
print(x)
y = np.arange(10, 20).reshape(1, 2, 5)
print(y)
z = tf.keras.layers.Dot(axes=(1, 2))([x, y])
z.shape
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# Loading data - select which cases to include in the training/validation set (commented out cases are held out)
cases = ['DUCT_1100',
         'DUCT_1150',
         'DUCT_1250',
         'DUCT_1300',
         'DUCT_1350',
         'DUCT_1400',
         'DUCT_1500',
         'DUCT_1600',
         'DUCT_1800',
         #'DUCT_2000',
         'DUCT_2205',
         'DUCT_2400',
         'DUCT_2600',
         'DUCT_2900',
         'DUCT_3200',
         #'DUCT_3500',
         'PHLL_case_0p5',
         'PHLL_case_0p8',
         'PHLL_case_1p0',
         #'PHLL_case_1p2',
         'PHLL_case_1p5',
         'BUMP_h20',
         'BUMP_h26',
         'BUMP_h31',
         #'BUMP_h38',
         'BUMP_h42',
         'CNDV_12600',
         'CNDV_20580',
         'CBFS_13700'
         ]

#Convenient functions for loading dataset
def loadCombinedArray(cases, field):
    # data = np.concatenate([
    #     np.load(f"/home/jupyter/ml-turbulence-dataset/{dataset}/{dataset}_{case}_{field}.npy", allow_pickle=True) 
    #     for case in cases])
    l = []
    for case in cases:
        path = f"/home/jupyter/momentum/ml-turbulence-dataset/{dataset}/{dataset}_{case}_{field}.npy"
        print(path)
        arr = np.load(path, allow_pickle=True)
        print(arr.shape)
        l.append(arr)
    return np.concatenate(l)

def loadLabels(cases, field):
    # data = np.concatenate([
    #     np.load(f"/home/jupyter/ml-turbulence-dataset/labels/{case}_{field}.npy", allow_pickle=True) 
    #     for case in cases])
    l = []
    for case in cases:
        path = f"/home/jupyter/momentum/ml-turbulence-dataset/labels/{case}_{field}.npy"
        print(path)
        arr = np.load(path, allow_pickle=True)
        print(arr.shape)
        l.append(arr)
    return np.concatenate(l)

# Select RANS model
dataset = 'kepsilon' 

print('Loading features and labels from the dataset: '+ dataset)

#Load the set of ten basis tensors (N,10,3,3), from Pope "A more general effective-viscosity hypothesis" (1975).
Tensors = loadCombinedArray(cases, 'Tensors')
print('Shape of basis tensor array: '+str(Tensors.shape))

#Load the 47 invariants (N,47) used by Wu et al. "Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework" (2018)
Invariants = loadCombinedArray(cases,'I1')
print('Shape of invariant features array: '+str(Invariants.shape))

#Load the additional scalars (N,5): 
#    q[:,0] = Ratio of excess rotation to strain rate,
#    q[:,1] = Wall-distance based Reynolds number, 
#    q[:,2] = Ratio of turbulent time scale to mean strain time scale
#    q[:,3] = Ratio of total Reynolds stress to 1/2 * normal Reynolds stress (TKE)
# Scalars = loadCombinedArray(cases,'q')
# print('Shape of scalar features array: '+str(Scalars.shape))

#Combine the invariants and scalars to form a feature set
# Features = np.column_stack((Invariants,Scalars))
Features = Invariants
print('Shape of combined features array: '+str(Features.shape))

#Optional: remove outliers based on the number of standard deviations away from the mean. 
#Note: must be careful, as there are naturally large variations in flow features. Even a 5*stdev critera removes many valid near-wall points.
def remove_outliers(Features):
    stdev = np.std(Features,axis=0)
    means = np.mean(Features,axis=0)
    ind_drop = np.empty(0)
    for i in range(len(Features[0,:])):
        ind_drop = np.concatenate((ind_drop,np.where((Features[:,i]>means[i]+5*stdev[i]) | (Features[:,i]<means[i]-5*stdev[i]) )[0]))
    return ind_drop.astype(int)

outlier_removal_switch = 0
if outlier_removal_switch == 1:
    outlier_index = remove_outliers(Features)
    print('Found '+str(len(outlier_index))+' outliers in the input feature set')
    Features = np.delete(Features,outlier_index,axis=0)
    Tensors = np.delete(Tensors,outlier_index,axis=0)
    Labels = np.delete(Labels,outlier_index,axis=0)

#Load the label set from DNS/LES:
Labels = loadLabels(cases,'b')
#If desired, reshape the 3x3 symmetric anisotropy tensor into a 1x6 vector
Labels = np.delete(Labels.reshape((len(Labels),9)),[3,6,7],axis=1)
print('Shape of DNS/LES labels array: '+str(Labels.shape))

# Split the datasets into training and validation
indices = np.arange(Features.shape[0])
input_shape = Features.shape[1]

x_train, x_val, y_train, y_val, ind_train, ind_val = train_test_split(Features, Labels, indices, test_size=0.2, random_state=10, shuffle=True)

basis_train = Tensors[ind_train]
basis_val = Tensors[ind_val]

scaler = StandardScaler()
x_train = scaler.fit_transform(x_train)
x_val = scaler.transform(x_val)

print(' ')
print('Training features:')
print(x_train.shape)
print('Training tensor basis:')
print(basis_train.shape)
print('Training labels:')
print(y_train.shape)
print(' ')
print('Validation features:')
print(x_val.shape)
print('Validation tensor basis:')
print(basis_val.shape)
print('Validation labels:')
print(y_val.shape)
print(' ')<jupyter_output>Loading features and labels from the dataset: kepsilon
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1100_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1150_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1250_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1300_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1350_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1400_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1500_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1600_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-turbulence-dataset/kepsilon/kepsilon_DUCT_1800_Tensors.npy
(9216, 10, 3, 3)
/home/jupyter/momentum/ml-t[...]<jupyter_text>Creating a tensor basis neural network (TBNN)The following model has 3 hidden layers, with 40 neurons each.<jupyter_code>import tensorflow as tf
from tensorflow import keras
keras.backend.clear_session()

#The model has two inputs, a set of input features with a learned mapping, and the tensor basis layer
input_layer = keras.layers.Input(shape=(input_shape),name ='input_layer')
input_tensor_basis = keras.layers.Input(shape=(10,3,3),name='Tensor_input_layer')

#Hidden layer definition
hidden1 = keras.layers.Dense(20,name='Hidden1', kernel_initializer="lecun_normal",kernel_regularizer=tf.keras.regularizers.l2(1E-3), activation = "selu")(input_layer)
hidden2 = keras.layers.Dense(20,name='Hidden2', kernel_initializer="lecun_normal",kernel_regularizer=tf.keras.regularizers.l2(1E-3), activation = "selu")(hidden1)
hidden3 = keras.layers.Dense(20,name='Hidden3', kernel_initializer="lecun_normal",kernel_regularizer=tf.keras.regularizers.l2(1E-3), activation = "selu")(hidden2)

#The layer of gn, which are coefficients for each of the ten Tn
gn = keras.layers.Dense(10,name='gn', kernel_initializer="lecun_normal",kernel_regularizer=tf.keras.regularizers.l2(1E-4), activation = "selu")(hidden3)
tf.print(gn.shape)

#Multiply the gn by Tn, with the output being the anisotropy tensor
shaped = keras.layers.Reshape((10,1,1),name='Shape_for_dot_product')(gn)
tf.print(shaped.shape)
tf.print(input_tensor_basis.shape)
merge = keras.layers.Dot(axes=1, name='Dot_product')([shaped,input_tensor_basis])
tf.print(f'merge shape {merge.shape}')

#Reshape the output anisotropy tensor, and trim out duplicate values (it is a symmetric matrix). The end result is a 6 component vector.
shaped_output = keras.layers.Reshape((9,1),name='Shaped_output')(merge)
tf.print(shaped_output.shape)
trimmed_output1 = keras.layers.Lambda(lambda x : x[:,0])(shaped_output)
trimmed_output2 = keras.layers.Lambda(lambda x : x[:,1])(shaped_output)
trimmed_output3 = keras.layers.Lambda(lambda x : x[:,2])(shaped_output)
trimmed_output4 = keras.layers.Lambda(lambda x : x[:,4])(shaped_output)
trimmed_output5 = keras.layers.Lambda(lambda x : x[:,5])(shaped_output)
trimmed_output6 = keras.layers.Lambda(lambda x : x[:,8])(shaped_output)
merged_output = tf.keras.layers.Concatenate()([trimmed_output1,trimmed_output2,trimmed_output3,trimmed_output4,trimmed_output5,trimmed_output6])
tf.print(merged_output.shape)

model=keras.Model(inputs=[input_layer, input_tensor_basis], outputs=[merged_output])


optimizer = tf.keras.optimizers.Nadam(learning_rate = 5E-4, clipnorm=1000)
model.compile(optimizer,loss='mse',metrics=['mae', 'mse'])
print(model.summary())<jupyter_output>TensorShape([None, 10])<jupyter_text>Train the tensor basis neural network<jupyter_code>reduce_lr =tf.keras.callbacks.ReduceLROnPlateau(
    monitor="loss",
    factor=0.3,
    patience=10,
    verbose=0,
    mode="auto",
    min_delta=0.0001,
    cooldown=0,
    min_lr=0,
)

early_stop = tf.keras.callbacks.EarlyStopping(
    monitor="val_loss",
    min_delta=0,
    patience=40,
    verbose=0,
    mode="auto",
    baseline=None,
    restore_best_weights=True,
)

history = model.fit([x_train,basis_train], y_train, 
                    batch_size=1000,
                    epochs=1000, 
                    validation_data = ([x_val,basis_val],y_val), 
                    verbose=1, 
                    callbacks=[
                            early_stop,
                            reduce_lr
                            ]
                   )
model.evaluate([x_train,basis_train],y_train)
model.evaluate([x_val,basis_val],y_val)
rand_ind = np.random.randint(0,len(x_train),5)
print('Picked 5 random indices: ' +str(rand_ind))

for i in range(len(rand_ind)):
    print('Index: '+str(rand_ind[i]))
    print('Label anisotropy values: ')
    print(y_train[rand_ind[i],:])
    print('Model prediction: ')
    print(model.predict([x_train[rand_ind[i],:].reshape(1,x_train.shape[1]),basis_train[rand_ind[i],:,:,:].reshape(1,10,3,3)]))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13171225937183384, 'fraction_numerical': 0.03564520585797182, 'mean_word_length': 4.522889114954221, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '533641', 'license': 'MIT', 'n_tokens_mistral': 4109, 'n_tokens_neox': 3807, 'n_words': 767, 'path': '_/example-use-of-ml-turbulence-dataset.ipynb', 'repo_name': 'alxyok/ransplacement', 'size': '10857'}