Feature multilayer perceptron and data-driven fluid model

.github/workflows/regression.yml

@@ -26,13 +26,13 @@ jobs:
         config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI, BaseOMP, ReverseOMP, ForwardOMP]
         include:
           - config_set: BaseMPI
-            flags: '-Denable-pywrapper=true -Denable-coolprop=true -Denable-tests=true --warnlevel=2 --werror'
+            flags: '-Denable-pywrapper=true -Denable-coolprop=true -Denable-mlpcpp=true -Denable-tests=true --warnlevel=2 --werror'
           - config_set: ReverseMPI
             flags: '-Denable-autodiff=true -Denable-normal=false -Denable-pywrapper=true -Denable-tests=true --warnlevel=3 --werror'
           - config_set: ForwardMPI
             flags: '-Denable-directdiff=true -Denable-normal=false -Denable-tests=true --warnlevel=3 --werror'
           - config_set: BaseNoMPI
-            flags: '-Denable-pywrapper=true -Denable-openblas=true -Dwith-mpi=disabled -Denable-tests=true --warnlevel=3 --werror'
+            flags: '-Denable-pywrapper=true -Denable-openblas=true -Dwith-mpi=disabled -Denable-mlpcpp=true -Denable-tests=true --warnlevel=3 --werror'
           - config_set: ReverseNoMPI
             flags: '-Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-pywrapper=true -Denable-tests=true --warnlevel=3 --werror'
           - config_set: ForwardNoMPI

.gitmodules

@@ -24,3 +24,6 @@
 [submodule "externals/mel"]
 	path = externals/mel
 	url = https://github.com/pcarruscag/MEL.git
+[submodule "subprojects/MLPCpp"]
+	path = subprojects/MLPCpp
+	url = https://github.com/EvertBunschoten/MLPCpp.git

Common/include/CConfig.hpp

@@ -519,7 +519,12 @@ class CConfig {
   Kind_TimeIntScheme_AdjTurb,   /*!< \brief Time integration for the adjoint turbulence model. */
   Kind_TimeIntScheme_Species,   /*!< \brief Time integration for the species model. */
   Kind_TimeIntScheme_Heat,      /*!< \brief Time integration for the wave equations. */
-  Kind_TimeStep_Heat;           /*!< \brief Time stepping method for the (fvm) heat equation. */
+  Kind_TimeStep_Heat,           /*!< \brief Time stepping method for the (fvm) heat equation. */
+  n_Datadriven_files;
+  ENUM_DATADRIVEN_METHOD Kind_DataDriven_Method;       /*!< \brief Method used for datset regression in data-driven fluid models. */
+
+  su2double DataDriven_Relaxation_Factor; /*!< \brief Relaxation factor for Newton solvers in data-driven fluid models. */
+
   STRUCT_TIME_INT Kind_TimeIntScheme_FEA;    /*!< \brief Time integration for the FEA equations. */
   STRUCT_SPACE_ITE Kind_SpaceIteScheme_FEA;  /*!< \brief Iterative scheme for nonlinear structural analysis. */
   unsigned short
@@ -800,7 +805,8 @@ class CConfig {
   SurfAdjCoeff_FileName,         /*!< \brief Output file with the adjoint variables on the surface. */
   SurfSens_FileName,             /*!< \brief Output file for the sensitivity on the surface (discrete adjoint). */
   VolSens_FileName,              /*!< \brief Output file for the sensitivity in the volume (discrete adjoint). */
-  ObjFunc_Hess_FileName;         /*!< \brief Hessian approximation obtained by the Sobolev smoothing solver. */
+  ObjFunc_Hess_FileName,         /*!< \brief Hessian approximation obtained by the Sobolev smoothing solver. */
+  *DataDriven_Method_FileNames;    /*!< \brief Dataset information for data-driven fluid models. */
 
   bool
   Wrt_Performance,           /*!< \brief Write the performance summary at the end of a calculation.  */
@@ -3844,6 +3850,29 @@ class CConfig {
    */
   unsigned short GetKind_FluidModel(void) const { return Kind_FluidModel; }
 
+  /*!
+   * \brief Datadriven method for EoS evaluation.
+   */
+  ENUM_DATADRIVEN_METHOD GetKind_DataDriven_Method(void) const { return Kind_DataDriven_Method; }
+
+  /*!
+   * \brief Get name of the input file for the data-driven fluid model interpolation method.
+   * \return Name of the input file for the interpolation method.
+   */
+  const string* GetDataDriven_FileNames(void) const { return DataDriven_Method_FileNames; }
+
+  /*!
+   * \brief Get number of listed look-up table or multi-layer perceptron input files.
+   * \return Number of listed data-driven method input files.
+   */
+  unsigned short GetNDataDriven_Files(void) const { return n_Datadriven_files; }
+
+  /*!
+   * \brief Get Newton solver relaxation factor for data-driven fluid models.
+   * \return Newton solver relaxation factor.
+   */
+  su2double GetRelaxation_DataDriven(void) const { return DataDriven_Relaxation_Factor; }
+
   /*!
    * \brief Returns the name of the fluid we are using in CoolProp.
    */

Common/include/option_structure.hpp

@@ -548,6 +548,7 @@ enum ENUM_FLUIDMODEL {
   FLUID_MIXTURE = 9,      /*!< \brief Species mixture model. */
   COOLPROP = 10,          /*!< \brief Thermodynamics library. */
   FLUID_FLAMELET = 11,    /*!< \brief lookup table (LUT) method for premixed flamelets. */
+  DATADRIVEN_FLUID = 12,           /*!< \brief multi-layer perceptron driven fluid model. */
 };
 static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("STANDARD_AIR", STANDARD_AIR)
@@ -561,6 +562,7 @@ static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("SU2_NONEQ", SU2_NONEQ)
   MakePair("FLUID_MIXTURE", FLUID_MIXTURE)
   MakePair("COOLPROP", COOLPROP)
+  MakePair("DATADRIVEN_FLUID", DATADRIVEN_FLUID)
   MakePair("FLUID_FLAMELET", FLUID_FLAMELET)
 };
 
@@ -590,6 +592,19 @@ MakePair("ARGON-SID",ARGON_SID)
 MakePair("ONESPECIES", ONESPECIES)
 };
 
+/*!
+* \brief Types of interpolation methods for data-driven fluid models.
+*/
+enum class ENUM_DATADRIVEN_METHOD {
+  LUT = 0,
+  MLP = 1
+};
+
+static const MapType<std::string, ENUM_DATADRIVEN_METHOD> DataDrivenMethod_Map = {
+  MakePair("LUT", ENUM_DATADRIVEN_METHOD::LUT)
+  MakePair("MLP", ENUM_DATADRIVEN_METHOD::MLP)
+};
+
 /*!
  * \brief types of coefficient transport model
  */

Common/src/CConfig.cpp

@@ -1169,6 +1169,14 @@ void CConfig::SetConfig_Options() {
   /*!\brief FLUID_NAME \n DESCRIPTION: Fluid name \n OPTIONS: see coolprop homepage \n DEFAULT: nitrogen \ingroup Config*/
   addStringOption("FLUID_NAME", FluidName, string("nitrogen"));
 
+  /*!\par CONFIG_CATEGORY: Data-driven fluid model parameters \ingroup Config*/
+  /*!\brief INTERPOLATION_METHOD \n DESCRIPTION: Interpolation method used to determine the thermodynamic state of the fluid. \n OPTIONS: See \link DataDrivenMethod_Map \endlink DEFAULT: MLP \ingroup Config*/
+  addEnumOption("INTERPOLATION_METHOD",Kind_DataDriven_Method, DataDrivenMethod_Map, ENUM_DATADRIVEN_METHOD::LUT);
+  /*!\brief FILENAME_INTERPOLATOR \n DESCRIPTION: Input file for the interpolation method. \n \ingroup Config*/
+  addStringListOption("FILENAMES_INTERPOLATOR", n_Datadriven_files, DataDriven_Method_FileNames);
+  /*!\brief DATADRIVEN_NEWTON_RELAXATION \n DESCRIPTION: Relaxation factor for Newton solvers in data-driven fluid model. \n \ingroup Config*/
+  addDoubleOption("DATADRIVEN_NEWTON_RELAXATION", DataDriven_Relaxation_Factor, 0.05);
+
   /*!\brief CONFINEMENT_PARAM \n DESCRIPTION: Input Confinement Parameter for Vorticity Confinement*/
   addDoubleOption("CONFINEMENT_PARAM", Confinement_Param, 0.0);
 
@@ -4967,6 +4975,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
   }
 
+  /*--- Data-driven fluid model is currently only supported for compressible flow problems. ---*/
+  if ((Kind_Solver == MAIN_SOLVER::INC_EULER || Kind_Solver == MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver == MAIN_SOLVER::INC_RANS) && (Kind_FluidModel == DATADRIVEN_FLUID)) {
+    SU2_MPI::Error("Data-driven fluid model can only be used for compressible flows.", CURRENT_FUNCTION);
+  }
+
   if ((Kind_Solver == MAIN_SOLVER::INC_EULER || Kind_Solver == MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver == MAIN_SOLVER::INC_RANS) && (Kind_FluidModel == INC_IDEAL_GAS_POLY)) {
     su2double sum = 0.0;
     for (unsigned short iVar = 0; iVar < N_POLY_COEFFS; iVar++) {

Docs/docmain.hpp

@@ -225,3 +225,9 @@
  * \brief Classes for explicit (done by the programmer) vectorization (SIMD) of computations.
  * \ingroup Toolboxes
  */
+
+/*!
+ * \defgroup Multi-Layer Perceptrons (MLP)
+ * \brief Data look up and interpolation via dense, feed-forward multi-layer perceptrons.
+ * \ingroup Toolboxes
+ */

SU2_CFD/include/fluid/CDataDrivenFluid.hpp

@@ -0,0 +1,220 @@
+/*!
+ * \file CDataDrivenFluid.hpp
+ * \brief Defines a template fluid model class using multilayer perceptrons
+ *  for theromodynamic state definition
+ * \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/>.
+ */
+
+#pragma once
+
+#include <vector>
+#include "../../../Common/include/containers/CLookUpTable.hpp"
+#if defined(HAVE_MLPCPP)
+#define MLP_CUSTOM_TYPE su2double
+#include "../../../subprojects/MLPCpp/include/CLookUp_ANN.hpp"
+#define USE_MLPCPP
+#endif
+#include "CFluidModel.hpp"
+
+/*!
+ * \class CDataDrivenFluid
+ * \brief Template class for fluid model definition using multi-layer perceptrons for
+ * fluid dynamic state definition.
+ * \author: E.C.Bunschoten.
+ */
+class CDataDrivenFluid final : public CFluidModel {
+ protected:
+  int rank{MASTER_NODE}; /*!< \brief MPI Rank. */
+  ENUM_DATADRIVEN_METHOD Kind_DataDriven_Method =
+      ENUM_DATADRIVEN_METHOD::LUT; /*!< \brief Interpolation method for data set evaluation. */
+
+  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. */
+      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. */
+
+  unsigned long MaxIter_Newton; /*!< \brief Maximum number of iterations for Newton solvers. */
+
+  su2double dsde_rho, /*!< \brief Entropy derivative w.r.t. density. */
+      dsdrho_e,       /*!< \brief Entropy derivative w.r.t. static energy. */
+      d2sde2,         /*!< \brief Entropy second derivative w.r.t. static energy. */
+      d2sdedrho,      /*!< \brief Entropy second derivative w.r.t. density and static energy. */
+      d2sdrho2;       /*!< \brief Entropy second derivative w.r.t. static density. */
+
+  su2double R_idealgas,     /*!< \brief Approximated ideal gas constant. */
+            Cp_idealgas,    /*!< \brief Approximated ideal gas specific heat at constant pressure. */
+            gamma_idealgas, /*!< \brief Approximated ideal gas specific heat ratio. */
+            Cv_idealgas,    /*!< \brief Approximated ideal gas specific heat at constant volume. */
+            P_middle,       /*!< \brief Pressure computed from the centre of the data set. */
+            T_middle;       /*!< \brief Temperature computed from the centre of the data set. */
+
+  su2double Enthalpy, /*!< \brief Fluid enthalpy value [J kg^-1] */
+      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<su2double*> outputs_rhoe; /*!< \brief Pointers to output variables. */
+
+  /*--- 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. */
+#endif
+  vector<su2double> MLP_inputs; /*!< \brief Inputs for the multi-layer perceptron look-up operation. */
+
+  CLookUpTable* lookup_table; /*!< \brief Look-up table regression object. */
+
+  unsigned long outside_dataset, /*!< \brief Density-energy combination lies outside data set. */
+      nIter_Newton;              /*!< \brief Number of Newton solver iterations. */
+
+  /*!
+   * \brief Map dataset variables to specific look-up operations.
+   */
+  void MapInputs_to_Outputs();
+
+  /*!
+   * \brief Evaluate dataset through multi-layer perceptron.
+   * \param[in] rho - Density value.
+   * \param[in] e - Static energy value.
+   * \return Query point lies outside MLP normalization range (0 is inside, 1 is outside).
+   */
+  unsigned long Predict_MLP(su2double rho, su2double e);
+
+  /*!
+   * \brief Evaluate dataset through look-up table.
+   * \param[in] rho - Density value.
+   * \param[in] e - Static energy value.
+   * \return Query point lies outside table data range (0 is inside, 1 is outside).
+   */
+  unsigned long Predict_LUT(su2double rho, su2double e);
+
+  /*!
+   * \brief Evaluate the data set.
+   * \param[in] rho - Density value.
+   * \param[in] e - Static energy value.
+   */
+  void Evaluate_Dataset(su2double rho, su2double e);
+
+  /*!
+   * \brief 2D Newton solver for computing the density and energy corresponding to Y1_target and Y2_target.
+   * \param[in] Y1_target - Target value for output quantity 1.
+   * \param[in] Y2_target - Target value for output quantity 2.
+   * \param[in] Y1 - Pointer to output quantity 1.
+   * \param[in] Y2 - Pointer to output quantity 2.
+   * \param[in] dY1drho - Pointer to the partial derivative of quantity 1 w.r.t. density at constant energy.
+   * \param[in] dY1de - Pointer to the partial derivative of quantity 1 w.r.t. energy at constant density.
+   * \param[in] dY2drho - Pointer to the partial derivative of quantity 2 w.r.t. density at constant energy.
+   * \param[in] dY2de - Pointer to the partial derivative of quantity 2 w.r.t. energy at constant density.
+   */
+  void Run_Newton_Solver(su2double Y1_target, su2double Y2_target, su2double* Y1, su2double* Y2, su2double* dY1drho,
+                         su2double* dY1de, su2double* dY2drho, su2double* dY2de);
+
+  /*!
+   * \brief 1D Newton solver for computing the density or energy corresponding to Y_target.
+   * \param[in] Y_target - Target quantity value.
+   * \param[in] Y - Pointer to output quantity.
+   * \param[in] X - Pointer to controlling variable (density or energy).
+   * \param[in] dYdX - Pointer to the partial derivative of target quantity w.r.t. controlling variable.
+   */
+  void Run_Newton_Solver(su2double Y_target, su2double* Y, su2double* X, su2double* dYdX);
+
+  void ComputeIdealGasQuantities();
+ public:
+  /*!
+   * \brief Constructor of the class.
+   */
+  CDataDrivenFluid(const CConfig* config, bool display = true);
+
+  ~CDataDrivenFluid();
+  /*!
+   * \brief Set the Dimensionless State using Density and Internal Energy.
+   * \param[in] rho - first thermodynamic variable (density).
+   * \param[in] e - second thermodynamic variable (static energy).
+   */
+  void SetTDState_rhoe(su2double rho, su2double e) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Pressure  and Temperature.
+   * \param[in] P - first thermodynamic variable (pressure).
+   * \param[in] T - second thermodynamic variable (temperature).
+   */
+  void SetTDState_PT(su2double P, su2double T) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Pressure and Density.
+   * \param[in] P - first thermodynamic variable (pressure).
+   * \param[in] rho - second thermodynamic variable (density).
+   */
+  void SetTDState_Prho(su2double P, su2double rho) override;
+
+  /*!
+   * \brief Set the Dimensionless Internal Energy using Pressure and Density.
+   * \param[in] P - first thermodynamic variable (pressure).
+   * \param[in] rho - second thermodynamic variable (density).
+   */
+  void SetEnergy_Prho(su2double P, su2double rho) override;
+
+  /*!
+   * \brief Set the Dimensionless Internal Energy using Pressure and Density.
+   * \param[in] rho - second thermodynamic variable (density).
+   */
+  void SetTDState_rhoT(su2double rho, su2double T) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Enthalpy and Entropy.
+   * \param[in] h - first thermodynamic variable (h).
+   * \param[in] s - second thermodynamic variable (s).
+   */
+  void SetTDState_hs(su2double h, su2double s) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Pressure and Entropy.
+   * \param[in] P - first thermodynamic variable (P).
+   * \param[in] s - second thermodynamic variable (s).
+   */
+  void SetTDState_Ps(su2double P, su2double s) override;
+
+  /*!
+   * \brief Get fluid model extrapolation instance.
+   * \return Query point lies outside fluid model data range.
+   */
+  unsigned long GetExtrapolation() override { return outside_dataset; }
+
+  /*!
+   * \brief Get number of Newton solver iterations.
+   * \return Newton solver iteration count at termination.
+   */
+  unsigned long GetnIter_Newton() override { return nIter_Newton; }
+};

SU2_CFD/include/fluid/CFluidModel.hpp

@@ -344,4 +344,16 @@ class CFluidModel {
    * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity.
    */
   void SetEddyViscosity(su2double val_Mu_Turb) { Mu_Turb = val_Mu_Turb; }
+
+  /*!
+   * \brief Get fluid model extrapolation instance
+   * \return Query point lies outside fluid model data range.
+   */
+  virtual unsigned long GetExtrapolation() { return 0; }
+
+  /*!
+   * \brief Get number of Newton solver iterations.
+   * \return Newton solver iteration count at termination.
+   */
+  virtual unsigned long GetnIter_Newton() { return 0; }
 };