Merge pull request #46 from su2code/develop

Develop

compressible_flow/NICFD_nozzle/DataDriven/Generate_Dataset.py

@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+
+## \file Generate_Dataset.py
+#  \brief Example python script for generating training data for 
+#   data-driven fluid model in SU2
+#  \author E.C.Bunschoten
+#  \version 7.5.1 "Blackbird"
+#
+# SU2 Project Website: https://su2code.github.io
+#
+# The SU2 Project is maintained by the SU2 Foundation
+# (http://su2foundation.org)
+#
+# Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md)
+#
+# SU2 is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# SU2 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+
+# make print(*args) function available in PY2.6+, does'nt work on PY < 2.6
+
+import CoolProp
+import numpy as np 
+from tqdm import tqdm
+import csv 
+
+# Name of the fluid in the CoolProp library.
+fluidName = 'Air'
+
+# Type of equation of state to be used by CoolProp.
+CP_eos = "HEOS"
+
+# Minimum and maximum dataset temperatures [K].
+T_min = 280
+T_max = 1000
+
+# Minimum and maximum dataset pressures [Pa].
+P_min = 5e4
+P_max = 2e6
+
+# Number of data points along each axis.
+Np_grid = 500
+
+# Fraction of data points to be used as training data for MLP training (0-1).
+f_train = 0.8
+
+# Fraction of data poins to be used as test data for MLP validation (0-1).
+f_test = 0.1
+
+
+# Prepare data grid
+T_range = np.linspace(T_min, T_max, Np_grid)
+P_range = np.linspace(P_min, P_max, Np_grid)
+
+T_grid, P_grid = np.meshgrid(T_range, P_range)
+
+T_dataset = T_grid.flatten()
+P_dataset = P_grid.flatten()
+
+density_dataset = np.zeros(np.shape(T_dataset))
+energy_dataset = np.zeros(np.shape(T_dataset))
+s_dataset = np.zeros(np.shape(T_dataset))
+dsde_dataset = np.zeros(np.shape(T_dataset))
+dsdrho_dataset = np.zeros(np.shape(T_dataset))
+d2sde2_dataset = np.zeros(np.shape(T_dataset))
+d2sdedrho_dataset = np.zeros(np.shape(T_dataset))
+d2sdrho2_dataset = np.zeros(np.shape(T_dataset))
+
+# Evaluate CoolProp on data grid.
+fluid = CoolProp.AbstractState(CP_eos, fluidName)
+idx_failed_below = []
+idx_failed_above = []
+print("Generating CoolProp data set...")
+for i in tqdm(range(len(T_dataset))):
+    try:
+        fluid.update(CoolProp.PT_INPUTS, P_dataset[i], T_dataset[i])
+
+        density_dataset[i] = fluid.rhomass()
+        energy_dataset[i] = fluid.umass()
+        s_dataset[i] = fluid.smass()
+        dsde_dataset[i] = fluid.first_partial_deriv(CoolProp.iSmass, CoolProp.iUmass, CoolProp.iDmass)
+        dsdrho_dataset[i] = fluid.first_partial_deriv(CoolProp.iSmass, CoolProp.iDmass, CoolProp.iUmass)
+        d2sde2_dataset[i] = fluid.second_partial_deriv(CoolProp.iSmass, CoolProp.iUmass, CoolProp.iDmass, CoolProp.iUmass, CoolProp.iDmass)
+        d2sdedrho_dataset[i] = fluid.second_partial_deriv(CoolProp.iSmass, CoolProp.iUmass, CoolProp.iDmass, CoolProp.iDmass, CoolProp.iUmass)
+        d2sdrho2_dataset[i] = fluid.second_partial_deriv(CoolProp.iSmass, CoolProp.iDmass, CoolProp.iUmass, CoolProp.iDmass, CoolProp.iUmass)
+    except:
+        idx_failed_below.append(i)
+        print("CoolProp failed at temperature "+str(T_dataset[i]) + ", pressure "+str(P_dataset[i]))
+print("Done!")
+
+# Collect all data arrays and fill in failed data points. 
+collected_data = np.vstack([density_dataset, 
+                                  energy_dataset, 
+                                  s_dataset, 
+                                  dsde_dataset, 
+                                  dsdrho_dataset, 
+                                  d2sde2_dataset, 
+                                  d2sdedrho_dataset, 
+                                  d2sdrho2_dataset]).T
+for i_failed in idx_failed_below:
+    collected_data[i_failed, :] = 0.5*(collected_data[i_failed+1, :] + collected_data[i_failed-1, :])
+
+# Shuffle data set and extract training, validation, and test data.
+np.random.shuffle(collected_data)
+np_train = int(f_train*len(density_dataset))
+np_val = int(f_test*len(density_dataset))
+np_test = len(density_dataset) - np_train - np_val
+
+train_data = collected_data[:np_train, :]
+dev_data = collected_data[np_train:(np_train+np_val), :]
+test_data = collected_data[(np_train+np_val):, :]
+
+# Write output files.
+with open(fluidName + "_dataset_full.csv", "w+") as fid:
+    fid.write("Density,Energy,s,dsde_rho,dsdrho_e,d2sde2,d2sdedrho,d2sdrho2\n")
+    csvWriter = csv.writer(fid,delimiter=',')
+    csvWriter.writerows(collected_data)
+
+with open(fluidName + "_dataset_train.csv", "w+") as fid:
+    fid.write("Density,Energy,s,dsde_rho,dsdrho_e,d2sde2,d2sdedrho,d2sdrho2\n")
+    csvWriter = csv.writer(fid,delimiter=',')
+    csvWriter.writerows(train_data)
+
+with open(fluidName + "_dataset_dev.csv", "w+") as fid:
+    fid.write("Density,Energy,s,dsde_rho,dsdrho_e,d2sde2,d2sdedrho,d2sdrho2\n")
+    csvWriter = csv.writer(fid,delimiter=',')
+    csvWriter.writerows(dev_data)
+
+with open(fluidName + "_dataset_test.csv", "w+") as fid:
+    fid.write("Density,Energy,s,dsde_rho,dsdrho_e,d2sde2,d2sdedrho,d2sdrho2\n")
+    csvWriter = csv.writer(fid,delimiter=',')
+    csvWriter.writerows(test_data)

compressible_flow/NICFD_nozzle/DataDriven/LUTWriter.m

@@ -0,0 +1,131 @@
+%!/usr/bin/env matlab
+
+%% \file LUTWriter.m
+%  \brief Example MATLAB script for generating a look-up table file 
+%   compatible with the CDataDriven_Fluid class in SU2.
+%  \author E.C.Bunschoten
+%  \version 7.5.1 "Blackbird"
+%
+% SU2 Project Website: https://su2code.github.io
+%
+% The SU2 Project is maintained by the SU2 Foundation
+% (http://su2foundation.org)
+%
+% Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md)
+%
+% SU2 is free software; you can redistribute it and/or
+% modify it under the terms of the GNU Lesser General Public
+% License as published by the Free Software Foundation; either
+% version 2.1 of the License, or (at your option) any later version.
+%
+% SU2 is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+% Lesser General Public License for more details.
+%
+% You should have received a copy of the GNU Lesser General Public
+% License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+
+% make print(*args) function available in PY2.6+, does'nt work on PY < 2.6
+
+% CoolProp input data file.
+input_datafile = "Air_dataset_full.csv";
+
+% Data point frequency (the larger, the coarser the table).
+data_freq = 2;
+
+% LUT output file name.
+output_LUTfile = "reftable.drg";
+
+%% provide import data file
+% space delimited
+data_ref = importdata(input_datafile);
+
+% Identify data entries
+rho = data_ref.data(1:data_freq:end, 1);
+e = data_ref.data(1:data_freq:end, 2);
+s = data_ref.data(1:data_freq:end, 3);
+ds_de = data_ref.data(1:data_freq:end, 4);
+ds_drho = data_ref.data(1:data_freq:end, 5);
+d2s_de2 = data_ref.data(1:data_freq:end, 6);
+d2s_dedrho = data_ref.data(1:data_freq:end, 7);
+d2s_drho2 = data_ref.data(1:data_freq:end, 8);
+
+rho_min = min(rho);
+rho_max = max(rho);
+e_min = min(e);
+e_max = max(e);
+
+% Normalize density and energy
+rho_norm = (rho - rho_min)/(rho_max - rho_min);
+e_norm = (e - e_min)/(e_max - e_min);
+
+%% Define table connectivity
+T = delaunayTriangulation(rho_norm, e_norm);
+
+data_LUT = [rho, e, s, ds_de, ds_drho, d2s_de2, d2s_dedrho, d2s_drho2];
+
+[~, boundNodes] = boundaryFacets(alphaShape(rho_norm, e_norm, 0.05));
+hullIDs = find(ismember([rho_norm, e_norm], boundNodes, "rows"));
+
+%% Write table data to output
+fid = fopen(output_LUTfile, 'w+');
+
+header = ['Dragon library' newline newline];
+
+header = [header '<Header>' newline];
+
+header = [header '[Version]' newline '1.0.1' newline newline];
+
+header = [header '[Number of points]' newline];
+header = [header sprintf('%3d',length(rho)) newline newline];
+
+header = [header '[Number of triangles]' newline];
+header = [header sprintf('%3d',length(T.ConnectivityList)) newline newline];
+
+header = [header '[Number of hull points]' newline];
+header = [header sprintf('%3d',length(hullIDs)) newline newline];
+
+header = [header '[Number of variables]' newline];
+header = [header sprintf('%3d',8) newline newline];
+
+header = [header '[Variable names]' newline];
+header = [header sprintf('1:Density\n2:Energy\n3:s\n4:dsde_rho\n5:dsdrho_e\n6:d2sde2\n7:d2sdedrho\n8:d2sdrho2\n')];
+
+header = [header newline '</Header>' newline newline];
+header = [header '<Data>'];
+
+fprintf(fid,'%s', header);
+printformat = '\n';
+for iTabVar=1:8
+    printformat = [printformat '%.14e\t'];
+end
+fprintf(fid,printformat,data_LUT');
+fprintf(fid,'%s', newline);
+fprintf(fid,'%s', '</Data>');
+
+fprintf(fid,'%s', newline);
+fprintf(fid,'%s', newline);
+fprintf(fid,'%s', '<Connectivity>');
+
+printformat = ['\n' '%5i\t' '%5i\t' '%5i\t'];
+
+fprintf(fid,printformat,T.ConnectivityList');
+
+fprintf(fid,'%s', newline);
+fprintf(fid,'%s', '</Connectivity>');
+
+%% print hull block
+fprintf(fid,'%s', newline);
+fprintf(fid,'%s', newline);
+fprintf(fid,'%s', '<Hull>');
+
+printformat = ['\n' '%5i\t'];
+
+fprintf(fid,printformat,hullIDs);
+
+fprintf(fid,'%s', newline);
+fprintf(fid,'%s', '</Hull>');
+
+%% close .dat file
+fclose(fid);