diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index c1bebd3928ff..cb7d6cca23ce 100644 --- a/.github/workflows/regression.yml +++ b/.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 diff --git a/SU2_CFD/include/fluid/CDataDrivenFluid.hpp b/SU2_CFD/include/fluid/CDataDrivenFluid.hpp index 45cd6f430625..47092a51f260 100644 --- a/SU2_CFD/include/fluid/CDataDrivenFluid.hpp +++ b/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 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 output_names_rhoe; /*!< \brief Output variable names listed in the data-driven method input file name. */ vector outputs_rhoe; /*!< \brief Pointers to output variables. */ vector> 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 MLP_inputs; /*!< \brief Inputs for the multi-layer perceptron look-up operation. */ diff --git a/SU2_CFD/src/fluid/CDataDrivenFluid.cpp b/SU2_CFD/src/fluid/CDataDrivenFluid.cpp index 64fa394f74f9..8b4f2c20859f 100644 --- a/SU2_CFD/src/fluid/CDataDrivenFluid.cpp +++ b/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: diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index e0d30ce4092e..d84a4b83a745 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/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& 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& 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; diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 622a1a6b068f..afc8ed5e2e52 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/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; diff --git a/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg b/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg index 93328224dc0b..974505d809bf 100644 --- a/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg +++ b/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 ---------------------------------------% diff --git a/UnitTests/Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp b/UnitTests/Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp index 2d255f8ca5cb..7f6dc38ad614 100644 --- a/UnitTests/Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp +++ b/UnitTests/Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp @@ -28,6 +28,7 @@ #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 @@ -35,43 +36,59 @@ #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 MLP_input_names, MLP_output_names; - std::vector MLP_inputs; - std::vector 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 input_names = {"x", "y"}, output_names = {"z"}; + std::vector output_refs = {&z_alt}; + std::vector 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 diff --git a/UnitTests/Common/toolboxes/multilayer_perceptron/MLP_Jacobian_tests.cpp b/UnitTests/Common/toolboxes/multilayer_perceptron/MLP_Jacobian_tests.cpp new file mode 100644 index 000000000000..b2277d7d1b5b --- /dev/null +++ b/UnitTests/Common/toolboxes/multilayer_perceptron/MLP_Jacobian_tests.cpp @@ -0,0 +1,119 @@ +/*! + * \file MLP_Jacobian_tests.cpp + * \brief Check if MLP is differentiable + * \author E.C.Bunschoten + * \version 8.4.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#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 + +#ifdef USE_MLPCPP +TEST_CASE("MLP Jacobian test", "[LookUpANN]") { + /* Create network with random weights */ + std::vector NN = {3, 10, 10, 1}; + MLPToolbox::CNeuralNetwork mlp = MLPToolbox::CNeuralNetwork(NN); + mlp.SetActivationFunction("gelu"); + mlp.RandomWeights(); + + /* Enable Jacobian and Hessian computation */ + mlp.CalcJacobian(true); + mlp.CalcHessian(true); + + /* Network input values */ + su2double val_in_1, val_in_2, val_in_3; + val_in_1 = 0.2; + val_in_2 = 0.3; + val_in_3 = 0.9; + + /* Calculate Jacobians with AD */ + AD::Reset(); + AD::StartRecording(); + AD::RegisterInput(val_in_1); + AD::RegisterInput(val_in_2); + AD::RegisterInput(val_in_3); + + mlp.SetInput(0, val_in_1); + mlp.SetInput(1, val_in_2); + mlp.SetInput(2, val_in_3); + mlp.Predict(); + + su2double val_out = mlp.GetOutput(0); + + AD::RegisterOutput(val_out); + AD::StopRecording(); + SU2_TYPE::SetDerivative(val_out, 1.0); + AD::ComputeAdjoint(); + const su2double jacobian_AD_1 = SU2_TYPE::GetDerivative(val_in_1); + const su2double jacobian_AD_2 = SU2_TYPE::GetDerivative(val_in_2); + const su2double jacobian_AD_3 = SU2_TYPE::GetDerivative(val_in_3); + + /* Retrieve analytical Jacobians */ + const su2double jacobian_ref_1 = mlp.GetJacobian(0, 0); + const su2double jacobian_ref_2 = mlp.GetJacobian(0, 1); + const su2double jacobian_ref_3 = mlp.GetJacobian(0, 2); + + /* Unit tests for Jacobians */ + CHECK(SU2_TYPE::GetValue(jacobian_AD_1) == Approx(SU2_TYPE::GetValue(jacobian_ref_1))); + CHECK(SU2_TYPE::GetValue(jacobian_AD_2) == Approx(SU2_TYPE::GetValue(jacobian_ref_2))); + CHECK(SU2_TYPE::GetValue(jacobian_AD_3) == Approx(SU2_TYPE::GetValue(jacobian_ref_3))); + + /* Calculate Hessians with AD */ + AD::Reset(); + AD::StartRecording(); + AD::RegisterInput(val_in_1); + AD::RegisterInput(val_in_2); + AD::RegisterInput(val_in_3); + + mlp.SetInput(0, val_in_1); + mlp.SetInput(1, val_in_2); + mlp.SetInput(2, val_in_3); + mlp.Predict(); + + su2double jacobian_analytical = mlp.GetJacobian(0, 0); + + AD::RegisterOutput(jacobian_analytical); + AD::StopRecording(); + SU2_TYPE::SetDerivative(jacobian_analytical, 1.0); + AD::ComputeAdjoint(); + const su2double hessian_AD_1 = SU2_TYPE::GetDerivative(val_in_1); + const su2double hessian_AD_2 = SU2_TYPE::GetDerivative(val_in_2); + const su2double hessian_AD_3 = SU2_TYPE::GetDerivative(val_in_3); + + /* Retrieve analytical Hessians */ + const su2double hessian_ref_1 = mlp.GetHessian(0, 0, 0); + const su2double hessian_ref_2 = mlp.GetHessian(0, 0, 1); + const su2double hessian_ref_3 = mlp.GetHessian(0, 0, 2); + + /* Unit tests on Hessians */ + CHECK(SU2_TYPE::GetValue(hessian_AD_1) == Approx(SU2_TYPE::GetValue(hessian_ref_1))); + CHECK(SU2_TYPE::GetValue(hessian_AD_2) == Approx(SU2_TYPE::GetValue(hessian_ref_2))); + CHECK(SU2_TYPE::GetValue(hessian_AD_3) == Approx(SU2_TYPE::GetValue(hessian_ref_3))); +} +#endif diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 3e64cfb9aa3d..40b54cf69717 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -19,10 +19,9 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp', 'SU2_CFD/windowing.cpp']) # Reverse-mode (algorithmic differentiation) tests: -su2_cfd_tests_ad = files(['Common/simple_ad_test.cpp']) -if get_option('enable-mlpcpp') - su2_cfd_tests_ad = su2_cfd_tests_ad + files(['SU2_CFD/fluid/CFluidModel_tests_AD.cpp']) -endif +su2_cfd_tests_ad = files(['Common/simple_ad_test.cpp', + 'SU2_CFD/fluid/CFluidModel_tests_AD.cpp', + 'Common/toolboxes/multilayer_perceptron/MLP_Jacobian_tests.cpp']) # Forward-mode (direct differentiation) tests: su2_cfd_tests_dd = files(['Common/simple_directdiff_test.cpp']) diff --git a/meson_scripts/init.py b/meson_scripts/init.py index 74e42b601776..8f54fd3b9594 100755 --- a/meson_scripts/init.py +++ b/meson_scripts/init.py @@ -74,7 +74,7 @@ def init_submodules( github_repo_mel = "https://github.com/pcarruscag/MEL" sha_version_fado = "ce7ee018e4e699af5028d69baa1939fea290e18a" github_repo_fado = "https://github.com/pcarruscag/FADO" - sha_version_mlpcpp = "ff57e0cf9e60127196d3f1be71e711d47ff646ef" + sha_version_mlpcpp = "02f2cb9dde791074858e11ac091f7c4df9c6af65" github_repo_mlpcpp = "https://github.com/EvertBunschoten/MLPCpp" sha_version_eigen = "d71c30c47858effcbd39967097a2d99ee48db464" github_repo_eigen = "https://gitlab.com/libeigen/eigen.git" diff --git a/subprojects/MLPCpp b/subprojects/MLPCpp index e19ca0cafb28..02f2cb9dde79 160000 --- a/subprojects/MLPCpp +++ b/subprojects/MLPCpp @@ -1 +1 @@ -Subproject commit e19ca0cafb28c4b7ba5b8cffef42883259b00dc0 +Subproject commit 02f2cb9dde791074858e11ac091f7c4df9c6af65