Update mlpcpp

.github/workflows/regression.yml

@@ -28,7 +28,7 @@ jobs:
           - config_set: BaseMPI
             flags: '-Denable-pywrapper=true -Denable-coolprop=true -Denable-mpp=true -Dinstall-mpp=true -Denable-mlpcpp=true -Denable-tests=true --warnlevel=2'
           - config_set: ReverseMPI
-            flags: '-Denable-autodiff=true -Denable-normal=false -Denable-pywrapper=true -Denable-tests=true --warnlevel=3 --werror'
+            flags: '-Denable-autodiff=true -Denable-normal=false -Denable-pywrapper=true -Denable-tests=true -Denable-mlpcpp=true --warnlevel=3 --werror'
           - config_set: ForwardMPI
             flags: '-Denable-directdiff=true -Denable-normal=false -Denable-tests=true -Denable-mlpcpp=true --warnlevel=3 --werror'
           - config_set: BaseNoMPI

SU2_CFD/include/fluid/CDataDrivenFluid.hpp

@@ -107,13 +107,12 @@ class CDataDrivenFluid final : public CFluidModel {
 
   string varname_rho,  /*!< \brief Controlling variable name for density. */
          varname_e;    /*!< \brief Controlling variable name for static energy. */
-
-  size_t idx_rho, /*!< \brief Interpolator index for density input. */
-      idx_e;      /*!< \brief Interpolator index for energy input. */
 
   su2double Newton_Relaxation, /*!< \brief Relaxation factor for Newton solvers. */
       rho_start,               /*!< \brief Starting value for the density in Newton solver processes. */
       e_start,                 /*!< \brief Starting value for the energy in Newton solver processes. */
+      rho_query,
+      e_query,
       Newton_Tolerance,        /*!< \brief Normalized tolerance for Newton solvers. */
       rho_min, rho_max,        /*!< \brief Minimum and maximum density values in data set. */
       e_min, e_max;            /*!< \brief Minimum and maximum energy values in data set. */
@@ -130,10 +129,7 @@ class CDataDrivenFluid final : public CFluidModel {
       dhdrho_e,       /*!< \brief Enthalpy derivative w.r.t. density. */
       dhde_rho;       /*!< \brief Enthalpy derivative w.r.t. static energy. */
 
-  vector<string> input_names_rhoe, /*!< \brief Data-driven method input variable names of the independent variables
-                                      (density, energy). */
-      output_names_rhoe; /*!< \brief Output variable names listed in the data-driven method input file name. */
-
+  vector<string> output_names_rhoe; /*!< \brief Output variable names listed in the data-driven method input file name. */
   vector<su2double*> outputs_rhoe; /*!< \brief Pointers to output variables. */
 
   vector<vector<su2double*>> dsdrhoe;           /*!< \brief Entropy Jacobian terms. */
@@ -144,8 +140,8 @@ class CDataDrivenFluid final : public CFluidModel {
   bool display_Newton_process{false};
   /*--- Class variables for the multi-layer perceptron method ---*/
 #ifdef USE_MLPCPP
-  MLPToolbox::CLookUp_ANN* lookup_mlp; /*!< \brief Multi-layer perceptron collection. */
-  MLPToolbox::CIOMap* iomap_rhoe;      /*!< \brief Input-output map. */
+  MLPToolbox::CLookUp_ANN *lookup_mlp; /*!< \brief Multi-layer perceptron collection. */
+  MLPToolbox::CIOMap iomap_rhoe;      /*!< \brief Input-output map. */
 #endif
   vector<su2double> MLP_inputs; /*!< \brief Inputs for the multi-layer perceptron look-up operation. */
 

SU2_CFD/src/fluid/CDataDrivenFluid.cpp

@@ -79,7 +79,6 @@ CDataDrivenFluid::~CDataDrivenFluid() {
   switch (Kind_DataDriven_Method) {
     case ENUM_DATADRIVEN_METHOD::MLP:
 #ifdef USE_MLPCPP
-      delete iomap_rhoe;
       delete lookup_mlp;
 #endif
       break;
@@ -92,82 +91,55 @@ CDataDrivenFluid::~CDataDrivenFluid() {
 }
 void CDataDrivenFluid::MapInputs_to_Outputs() {
   /*--- Inputs of the data-driven method are density and internal energy. ---*/
-  input_names_rhoe.resize(2);
-  idx_rho = 0;
-  idx_e = 1;
-  input_names_rhoe[idx_rho] = varname_rho;
-  input_names_rhoe[idx_e] = varname_e;
 
   /*--- Required outputs for the interpolation method are entropy and its partial derivatives with respect to energy and
    * density. ---*/
-  size_t n_outputs, idx_s,idx_dsde_rho = 1, idx_dsdrho_e = 2, idx_d2sde2 = 3, idx_d2sdedrho = 4, idx_d2sdrho2 = 5;
-  if (use_MLP_derivatives) {
-    n_outputs = 1;
-    idx_s = 0;
-
-    outputs_rhoe.resize(n_outputs);
-    output_names_rhoe.resize(n_outputs);
-    output_names_rhoe[idx_s] = "s";
-    outputs_rhoe[idx_s] = &Entropy;
-
-    dsdrhoe.resize(n_outputs);
-    d2sdrhoe2.resize(n_outputs);
-    dsdrhoe[0].resize(2);
-    dsdrhoe[0][idx_rho] = &dsdrho_e;
-    dsdrhoe[0][idx_e] = &dsde_rho;
-
-    d2sdrhoe2[0].resize(2);
-    d2sdrhoe2[0][idx_rho].resize(2);
-    d2sdrhoe2[0][idx_e].resize(2);
-    d2sdrhoe2[0][idx_rho][idx_rho] = &d2sdrho2;
-    d2sdrhoe2[0][idx_rho][idx_e] = &d2sdedrho;
-    d2sdrhoe2[0][idx_e][idx_rho] = &d2sdedrho;
-    d2sdrhoe2[0][idx_e][idx_e] = &d2sde2;
-  } else {
-    n_outputs = 6;
-    idx_s = 0;
-    idx_dsde_rho = 1, idx_dsdrho_e = 2, idx_d2sde2 = 3, idx_d2sdedrho = 4, idx_d2sdrho2 = 5;
-
-    outputs_rhoe.resize(n_outputs);
-    output_names_rhoe.resize(n_outputs);
-    output_names_rhoe[idx_s] = "s";
-    outputs_rhoe[idx_s] = &Entropy;
-    output_names_rhoe[idx_dsde_rho] = "dsde_rho";
-    outputs_rhoe[idx_dsde_rho] = &dsde_rho;
-    output_names_rhoe[idx_dsdrho_e] = "dsdrho_e";
-    outputs_rhoe[idx_dsdrho_e] = &dsdrho_e;
-    output_names_rhoe[idx_d2sde2] = "d2sde2";
-    outputs_rhoe[idx_d2sde2] = &d2sde2;
-    output_names_rhoe[idx_d2sdedrho] = "d2sdedrho";
-    outputs_rhoe[idx_d2sdedrho] = &d2sdedrho;
-    output_names_rhoe[idx_d2sdrho2] = "d2sdrho2";
-    outputs_rhoe[idx_d2sdrho2] = &d2sdrho2;
-  }
-
+  const string name_entropy{"s"}, name_dsdrho{"dsdrho_e"}, name_dsde{"dsde"}, name_d2sdrho2{"d2sdrho2"}, name_d2sdedrho{"d2sdedrho"}, name_d2sde2{"d2sde2"};
+
 
-  /*--- Further preprocessing of input and output variables. ---*/
   if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::MLP) {
-/*--- Map MLP inputs to outputs. ---*/
-#ifdef USE_MLPCPP
-    iomap_rhoe = new MLPToolbox::CIOMap(input_names_rhoe, output_names_rhoe);
-    lookup_mlp->PairVariableswithMLPs(*iomap_rhoe);
-    MLP_inputs.resize(2);
-#endif
+    #ifdef USE_MLPCPP
+    iomap_rhoe = MLPToolbox::CIOMap();
+    iomap_rhoe.AddQueryInput(varname_rho, &rho_query);
+    iomap_rhoe.AddQueryInput(varname_e, &e_query);
+    iomap_rhoe.AddQueryOutput(name_entropy, &Entropy);
+
+    if (use_MLP_derivatives) {
+      iomap_rhoe.AddQueryJacobian(name_entropy, varname_rho, &dsdrho_e);
+      iomap_rhoe.AddQueryJacobian(name_entropy, varname_e, &dsde_rho);
+      iomap_rhoe.AddQueryHessian(name_entropy, varname_rho, varname_rho, &d2sdrho2);
+      iomap_rhoe.AddQueryHessian(name_entropy, varname_e, varname_rho, &d2sdedrho);
+      iomap_rhoe.AddQueryHessian(name_entropy, varname_e, varname_e, &d2sde2);
+    } else {
+      iomap_rhoe.AddQueryOutput(name_dsdrho, &dsdrho_e);
+      iomap_rhoe.AddQueryOutput(name_dsde, &dsde_rho);
+      iomap_rhoe.AddQueryOutput(name_d2sdrho2, &d2sdrho2);
+      iomap_rhoe.AddQueryOutput(name_d2sde2, &d2sde2);
+      iomap_rhoe.AddQueryOutput(name_d2sdedrho, &d2sdedrho);
+    }
+
+    lookup_mlp->PairVariableswithMLPs(iomap_rhoe);
+    #endif
   } else {
-    /*--- Retrieve column indices of LUT output variables ---*/
-    LUT_idx_s = lookup_table->GetIndexOfVar(output_names_rhoe[idx_s]);
-    LUT_idx_dsdrho_e = lookup_table->GetIndexOfVar(output_names_rhoe[idx_dsdrho_e]);
-    LUT_idx_dsde_rho = lookup_table->GetIndexOfVar(output_names_rhoe[idx_dsde_rho]);
-    LUT_idx_d2sde2 = lookup_table->GetIndexOfVar(output_names_rhoe[idx_d2sde2]);
-    LUT_idx_d2sdedrho= lookup_table->GetIndexOfVar(output_names_rhoe[idx_d2sdedrho]);
-    LUT_idx_d2sdrho2 = lookup_table->GetIndexOfVar(output_names_rhoe[idx_d2sdrho2]);
+    LUT_idx_s = lookup_table->GetIndexOfVar(name_entropy);
+    LUT_idx_dsdrho_e = lookup_table->GetIndexOfVar(name_dsdrho);
+    LUT_idx_dsde_rho = lookup_table->GetIndexOfVar(name_dsde);
+    LUT_idx_d2sde2 = lookup_table->GetIndexOfVar(name_d2sde2);
+    LUT_idx_d2sdedrho= lookup_table->GetIndexOfVar(name_d2sdedrho);
+    LUT_idx_d2sdrho2 = lookup_table->GetIndexOfVar(name_d2sdrho2);
 
     LUT_lookup_indices.push_back(LUT_idx_s);
+    outputs_rhoe.push_back(&Entropy);
     LUT_lookup_indices.push_back(LUT_idx_dsde_rho);
+    outputs_rhoe.push_back(&dsde_rho);
     LUT_lookup_indices.push_back(LUT_idx_dsdrho_e);
+    outputs_rhoe.push_back(&dsdrho_e);
     LUT_lookup_indices.push_back(LUT_idx_d2sde2);
+    outputs_rhoe.push_back(&d2sde2);
     LUT_lookup_indices.push_back(LUT_idx_d2sdedrho);
+    outputs_rhoe.push_back(&d2sdedrho);
     LUT_lookup_indices.push_back(LUT_idx_d2sdrho2);
+    outputs_rhoe.push_back(&d2sdrho2);
   }
 }
 
@@ -297,13 +269,9 @@ unsigned long CDataDrivenFluid::Predict_MLP(su2double rho, su2double e) {
   unsigned long exit_code = 0;
 /*--- Evaluate MLP collection for the given values for density and energy. ---*/
 #ifdef USE_MLPCPP
-  MLP_inputs[idx_rho] = rho;
-  MLP_inputs[idx_e] = e;
-  if (use_MLP_derivatives){
-    exit_code = lookup_mlp->PredictANN(iomap_rhoe, MLP_inputs, outputs_rhoe, &dsdrhoe, &d2sdrhoe2);
-  } else {
-    exit_code = lookup_mlp->PredictANN(iomap_rhoe, MLP_inputs, outputs_rhoe);
-  }
+  rho_query = rho;
+  e_query = e;
+  if(!lookup_mlp->Predict(iomap_rhoe)) exit_code=1;
 #endif
 
   return exit_code;
@@ -435,10 +403,10 @@ void CDataDrivenFluid::ComputeIdealGasQuantities() {
     break;
   case ENUM_DATADRIVEN_METHOD::MLP:
 #ifdef USE_MLPCPP
-    rho_min = lookup_mlp->GetInputNorm(iomap_rhoe, idx_rho).first;
-    e_min = lookup_mlp->GetInputNorm(iomap_rhoe, idx_e).first;
-    rho_max = lookup_mlp->GetInputNorm(iomap_rhoe, idx_rho).second;
-    e_max = lookup_mlp->GetInputNorm(iomap_rhoe, idx_e).second;
+    rho_min = iomap_rhoe.GetInputNorm(varname_rho).first;
+    e_min = iomap_rhoe.GetInputNorm(varname_e).first;
+    rho_max = iomap_rhoe.GetInputNorm(varname_rho).second;
+    e_max = iomap_rhoe.GetInputNorm(varname_e).second;
 #endif
     break;
   default:

SU2_CFD/src/fluid/CFluidFlamelet.cpp

@@ -224,23 +224,6 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) {
   val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC] = beta_mixfrac;
 
   preferential_diffusion = flamelet_options.preferential_diffusion;
-  switch (Kind_DataDriven_Method) {
-    case ENUM_DATADRIVEN_METHOD::LUT:
-      preferential_diffusion = look_up_table->CheckForVariables(varnames_PD);
-      break;
-    case ENUM_DATADRIVEN_METHOD::MLP:
-#ifdef USE_MLPCPP
-      n_betas = 0;
-      for (auto iMLP = 0u; iMLP < datadriven_fluid_options.n_filenames; iMLP++) {
-        auto outputMap = lookup_mlp->FindVariableIndices(iMLP, varnames_PD, false);
-        n_betas += outputMap.size();
-      }
-      preferential_diffusion = (n_betas == varnames_PD.size());
-#endif
-      break;
-    default:
-      break;
-  }
 
   if (!preferential_diffusion && flamelet_options.preferential_diffusion)
     SU2_MPI::Error("Preferential diffusion scalars not included in flamelet manifold.", CURRENT_FUNCTION);
@@ -252,7 +235,8 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) {
     iomap_LookUp = new MLPToolbox::CIOMap(controlling_variable_names, varnames_LookUp);
     lookup_mlp->PairVariableswithMLPs(*iomap_TD);
     lookup_mlp->PairVariableswithMLPs(*iomap_Sources);
-    lookup_mlp->PairVariableswithMLPs(*iomap_LookUp);
+    if (n_lookups > 1)
+      lookup_mlp->PairVariableswithMLPs(*iomap_LookUp);
     if (preferential_diffusion) {
       iomap_PD = new MLPToolbox::CIOMap(controlling_variable_names, varnames_PD);
       lookup_mlp->PairVariableswithMLPs(*iomap_PD);
@@ -329,7 +313,7 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector<su2double>& input_sca
 
 
   /*--- Add all quantities and their names to the look up vectors. ---*/
-  bool inside;
+  bool inside{true};
   switch (Kind_DataDriven_Method) {
     case ENUM_DATADRIVEN_METHOD::LUT:
       if (output_refs.size() != LUT_idx.size())
@@ -339,19 +323,20 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector<su2double>& input_sca
       } else {
         inside = look_up_table->LookUp_XY(LUT_idx, output_refs, val_prog, val_enth);
       }
-      if (inside) extrapolation = 0;
-      else extrapolation = 1;
+
       break;
     case ENUM_DATADRIVEN_METHOD::MLP:
       refs_vars.resize(output_refs.size());
       for (auto iVar = 0u; iVar < output_refs.size(); iVar++) refs_vars[iVar] = &output_refs[iVar];
 #ifdef USE_MLPCPP
-      extrapolation = lookup_mlp->PredictANN(iomap_Current, input_scalar, refs_vars);
+      inside=lookup_mlp->Predict(*iomap_Current, input_scalar, refs_vars);
 #endif
       break;
     default:
       break;
   }
+  if (inside) extrapolation = 0;
+      else extrapolation = 1;
   for (auto iVar = 0u; iVar < output_refs.size(); iVar++) AD::SetPreaccOut(output_refs[iVar]);
   AD::EndPreacc();
   return extrapolation;

SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp

@@ -213,7 +213,8 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver**
 
     for (unsigned long i_mesh = 0; i_mesh <= config->GetnMGLevels(); i_mesh++) {
       fluid_model_local = solver_container[i_mesh][FLOW_SOL]->GetFluidModel();
-      prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init);
+      if (flame_front_ignition) prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init);
+
       for (auto iVar = 0u; iVar < nVar; iVar++) scalar_init[iVar] = config->GetSpecies_Init()[iVar];
 
       /*--- Set enthalpy based on initial temperature and scalars. ---*/
@@ -797,7 +798,7 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo
   bool outside = false;
   while (!outside) {
     fluid_model->SetTDState_T(300, scalars);
-    if (fluid_model->GetExtrapolation() == 1) outside = true;
+    if (fluid_model->GetExtrapolation() == 1 || (fluid_model->GetTemperature()>1000.)) outside = true;
     scalars[I_PROGVAR] += delta;
   }
   su2double pv_burnt = scalars[I_PROGVAR] - delta;

TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg

@@ -31,7 +31,7 @@ CONDUCTIVITY_MODEL= FLAMELET
 DIFFUSIVITY_MODEL= FLAMELET
 KIND_SCALAR_MODEL= FLAMELET
 INTERPOLATION_METHOD= MLP
-FILENAMES_INTERPOLATOR= (MLP_TD.mlp, MLP_PD.mlp, MLP_PPV.mlp, MLP_null.mlp)
+FILENAMES_INTERPOLATOR= (MLP_TD.mlp, MLP_PD.mlp, MLP_PPV.mlp)
 PREFERENTIAL_DIFFUSION= YES
 
 % -------------------- SCALAR TRANSPORT ---------------------------------------%

UnitTests/Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp

@@ -28,50 +28,67 @@
 #include "catch.hpp"
 #include "../../../../Common/include/CConfig.hpp"
 #if defined(HAVE_MLPCPP)
+#define MLP_CUSTOM_TYPE su2double
 #include "../../../../subprojects/MLPCpp/include/CLookUp_ANN.hpp"
 #define USE_MLPCPP
 #endif
 #include <vector>
 
 #ifdef USE_MLPCPP
 TEST_CASE("LookUp ANN test", "[LookUpANN]") {
-  std::string MLP_input_files[] = {"src/SU2/UnitTests/Common/toolboxes/multilayer_perceptron/simple_mlp.mlp"};
-  unsigned short n_MLPs = 1;
-  MLPToolbox::CLookUp_ANN ANN(n_MLPs, MLP_input_files);
-  std::vector<std::string> MLP_input_names, MLP_output_names;
-  std::vector<double> MLP_inputs;
-  std::vector<double*> MLP_outputs;
-  su2double x, y, z;
+  MLPToolbox::CLookUp_ANN ANN;
+  MLPToolbox::CNeuralNetwork mlp =
+      MLPToolbox::CNeuralNetwork("src/SU2/UnitTests/Common/toolboxes/multilayer_perceptron/simple_mlp.mlp");
 
-  /*--- Define MLP inputs and outputs ---*/
-  MLP_input_names.resize(2);
-  MLP_input_names[0] = "x";
-  MLP_input_names[1] = "y";
-  MLP_inputs.resize(2);
+  ANN.AddNetwork(&mlp);
+  su2double x, y, z, z_alt;
 
-  MLP_outputs.resize(1);
-  MLP_output_names.resize(1);
-  MLP_output_names[0] = "z";
-  MLP_outputs[0] = &z;
+  /*--- Create a query where the value of z is calculated from x and y ---*/
+  MLPToolbox::CIOMap iomap_ref;
+  iomap_ref.AddQueryInput("x", &x);
+  iomap_ref.AddQueryInput("y", &y);
+  iomap_ref.AddQueryOutput("z", &z);
 
-  /*--- Generate input-output map ---*/
-  MLPToolbox::CIOMap iomap(MLP_input_names, MLP_output_names);
-  ANN.PairVariableswithMLPs(iomap);
-  /*--- MLP evaluation on point in the middle of the training data range ---*/
+  ANN.PairVariableswithMLPs(iomap_ref);
+
+  MLPToolbox::CIOMap iomap_vec;
+  std::vector<std::string> input_names = {"x", "y"}, output_names = {"z"};
+  std::vector<su2double*> output_refs = {&z_alt};
+  std::vector<su2double> mlp_inputs;
+
+  iomap_vec.SetQueryInput(input_names);
+  iomap_vec.SetQueryOutput(output_names);
+  ANN.PairVariableswithMLPs(iomap_vec);
+
+  /*--- MLP evaluation on two points in the middle of the training data range ---*/
   x = 1.0;
   y = -0.5;
 
-  MLP_inputs[0] = x;
-  MLP_inputs[1] = y;
-  ANN.PredictANN(&iomap, MLP_inputs, MLP_outputs);
+  bool inside = ANN.Predict(iomap_ref);
   CHECK(z == Approx(0.344829));
+  CHECK(inside);
+
+  mlp_inputs.resize(2);
+  mlp_inputs[0] = x;
+  mlp_inputs[1] = y;
+  ANN.Predict(iomap_vec, mlp_inputs, output_refs);
+  CHECK(z == z_alt);
+
+  x = 0.5;
+  y = -0.23;
+  inside = ANN.Predict(iomap_ref);
+
+  CHECK(z == Approx(0.224986));
+  CHECK(inside);
+
+  mlp_inputs[0] = x;
+  mlp_inputs[1] = y;
+  ANN.Predict(iomap_vec, mlp_inputs, output_refs);
 
-  /*--- MLP evaluation on point outside the training data range ---*/
-  x = 3.0;
-  y = -10;
-  MLP_inputs[0] = x;
-  MLP_inputs[1] = y;
-  ANN.PredictANN(&iomap, MLP_inputs, MLP_outputs);
-  CHECK(z == Approx(0.012737));
+  /*--- */
+  x = 10.0;
+  y = -20.0;
+  inside = ANN.Predict(iomap_ref);
+  CHECK(!inside);
 }
 #endif