diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 31aaed6f6930..dde6c94cefff 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -211,7 +211,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c fix_flame_initialization -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -260,7 +260,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests" + args: -b ${{github.ref}} -t develop -c fix_flame_initialization -s ${{matrix.testscript}} -a "--tapetests" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -308,7 +308,7 @@ jobs: PMIX_MCA_gds: hash with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan" + args: -b ${{github.ref}} -t develop -c fix_flame_initialization -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:260405-0054 with: @@ -353,7 +353,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan" + args: -b ${{github.ref}} -t develop -c fix_flame_initialization-s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index b3d2478cf5f2..477449031912 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1461,6 +1461,8 @@ struct FluidFlamelet_ParsedOptions { su2double* spark_reaction_rates; /*!< \brief Source terms for flamelet spark ignition option. */ unsigned short nspark; /*!< \brief Number of source terms for spark initialization. */ bool preferential_diffusion = false; /*!< \brief Preferential diffusion physics for flamelet solver.*/ + su2double Flame_T_ignition = 5000; /*!< \brief Ignition temperature for the flame, used for initialization. */ + }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 2ba453b197b0..c2b2f8ff9cbd 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1419,6 +1419,9 @@ void CConfig::SetConfig_Options() { /*!\brief SPARK_REACTION_RATES \n DESCRIPTION: Net source term values applied to species within spark area during spark ignition. \ingroup Config*/ addDoubleListOption("SPARK_REACTION_RATES", flamelet_ParsedOptions.nspark, flamelet_ParsedOptions.spark_reaction_rates); + /*!\brief FLAME_INIT_IGNITION \n DESCRIPTION: Ignition temperature for the flame initialization \ingroup Config*/ + addDoubleOption("FLAME_INIT_IGNITION", flamelet_ParsedOptions.Flame_T_ignition, 5000.0); + /*--- Options related to mass diffusivity and thereby the species solver. ---*/ /*!\brief DIFFUSIVITY_MODEL\n DESCRIPTION: mass diffusivity model \n DEFAULT constant disffusivity \ingroup Config*/ @@ -5700,6 +5703,26 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("Number of initial species incompatible with number of controlling variables and user scalars.", CURRENT_FUNCTION); /*--- We can have additional user defined transported scalars ---*/ flamelet_ParsedOptions.n_scalars = flamelet_ParsedOptions.n_control_vars + flamelet_ParsedOptions.n_user_scalars; + + /*--- Check that spark ignition has required parameters defined ---*/ + if (flamelet_ParsedOptions.ignition_method == FLAMELET_INIT_TYPE::SPARK) { + /*--- Check if SPARK_INIT was explicitly set in config file ---*/ + if (all_options.find("SPARK_INIT") != all_options.end()) { + SU2_MPI::Error("FLAME_INIT_METHOD= SPARK requires SPARK_INIT to be defined in the config file.", CURRENT_FUNCTION); + } + /*--- Check if SPARK_REACTION_RATES was explicitly set in config file ---*/ + if (all_options.find("SPARK_REACTION_RATES") != all_options.end()) { + SU2_MPI::Error("FLAME_INIT_METHOD= SPARK requires SPARK_REACTION_RATES to be defined in the config file.", CURRENT_FUNCTION); + } + if (flamelet_ParsedOptions.nspark < flamelet_ParsedOptions.n_scalars) { + SU2_MPI::Error("SPARK_REACTION_RATES must have at least " + to_string(flamelet_ParsedOptions.n_scalars) + + " values (one for each scalar variable), but only " + to_string(flamelet_ParsedOptions.nspark) + " were provided.", CURRENT_FUNCTION); + } + } + /*--- Check if flame ignition temperature is valid ---*/ + if (flamelet_ParsedOptions.Flame_T_ignition <= Inc_Temperature_Init) { + SU2_MPI::Error("Flame ignition temperature must be higher than the initial temperature of the flow field.", CURRENT_FUNCTION); + } } if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) { diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index df61813a3e27..5801cd6d8e3c 100644 --- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp @@ -70,9 +70,10 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { * \brief Find maximum progress variable value within the manifold for the current solution. * \param[in] fluid_model - pointer to flamelet fluid model. * \param[in] scalars - local scalar solution. + * \param[in] T_ignition - ignition temperature - at this temperature, the progress variable is considered burnt. * \return - maximum progress variable value within manifold bounds. */ - su2double GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalars); + su2double GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalars, const su2double T_ignition); /*! * \brief Retrieve scalar source terms from manifold. diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index de8d5e618531..297c1b1540b0 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -103,9 +103,13 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver su2double dist_from_center = 0, spark_radius = flamelet_config_options.spark_init[3]; dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), flamelet_config_options.spark_init.data()); - if (dist_from_center < pow(spark_radius,2)) { - for (auto iVar = 0u; iVar < nVar; iVar++) - nodes->SetScalarSource(i_point, iVar, nodes->GetScalarSources(i_point)[iVar] + flamelet_config_options.spark_reaction_rates[iVar]); + su2double T_local = flowNodes->GetTemperature(i_point); + if (dist_from_center < pow(spark_radius,2) && T_local < flamelet_config_options.Flame_T_ignition) { + /*--- Add spark reaction rates to the sources that were just set by SetScalarSources ---*/ + const su2double* current_sources = nodes->GetScalarSources(i_point); + for (auto iVar = 0u; iVar < nVar; iVar++) { + nodes->SetScalarSource(i_point, iVar, current_sources[iVar] + flamelet_config_options.spark_reaction_rates[iVar]); + } } } @@ -184,7 +188,7 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** su2double enth_inlet = config->GetSpecies_Init()[I_ENTH]; su2double prog_burnt = 0, prog_unburnt, point_loc; - su2double scalar_init[MAXNVAR]; + su2double scalar_init[MAXNVAR]= {0.0}; if (rank == MASTER_NODE) { cout << "initial condition: T = " << temp_inlet << endl; @@ -197,10 +201,11 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** cout << "Ignition with a straight flame front" << endl; break; case FLAMELET_INIT_TYPE::SPARK: - cout << "Ignition with an artificial spark" << endl; + cout << "Ignition with an artificial spark at iteration "<< flamelet_config_options.spark_init[4] + << " for a duration of " << flamelet_config_options.spark_init[5] << " iterations." << endl; break; case FLAMELET_INIT_TYPE::NONE: - cout << "No solution ignition (cold flow)" << endl; + cout << "No solution ignition (cold flow or restart)" << endl; break; default: break; @@ -216,7 +221,6 @@ 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(); - 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]; @@ -224,6 +228,8 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** n_not_iterated_local += GetEnthFromTemp(fluid_model_local, temp_inlet, config->GetSpecies_Init(), &enth_inlet); scalar_init[I_ENTH] = enth_inlet; + if (flame_front_ignition) prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init, flamelet_config_options.Flame_T_ignition); + prog_unburnt = config->GetSpecies_Init()[I_PROGVAR]; SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long i_point = 0; i_point < nPoint; i_point++) { @@ -805,16 +811,22 @@ unsigned long CSpeciesFlameletSolver::GetEnthFromTemp(CFluidModel* fluid_model, return exit_code; } -su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalar_solution) { +su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_model, const su2double* scalar_solution, const su2double T_ignition) { SU2_ZONE_SCOPED su2double scalars[MAXNVAR], delta = 1e-3; for (auto iVar = 0u; iVar < nVar; iVar++) scalars[iVar] = scalar_solution[iVar]; bool outside = false; + scalars[I_PROGVAR] += delta; while (!outside) { - fluid_model->SetTDState_T(300, scalars); - if (fluid_model->GetExtrapolation() == 1 || (fluid_model->GetTemperature()>1000.)) outside = true; + /*--- Note that 300.0 is a dummy temperature here and not used. ---*/ + fluid_model->SetTDState_T(300.0, scalars); + if ((fluid_model->GetExtrapolation() == 1) || fluid_model->GetTemperature() > T_ignition) outside = true; scalars[I_PROGVAR] += delta; } su2double pv_burnt = scalars[I_PROGVAR] - delta; + if (rank == MASTER_NODE) { + cout << "Burnt progress variable determined from flamelet table: " << pv_burnt << endl; + cout << "Burnt temperature from flamelet table: " << fluid_model->GetTemperature() << endl; + } return pv_burnt; } diff --git a/TestCases/flamelet/08_flame_init_method/flamelet_initialization.py b/TestCases/flamelet/08_flame_init_method/flamelet_initialization.py new file mode 100644 index 000000000000..514a9f140127 --- /dev/null +++ b/TestCases/flamelet/08_flame_init_method/flamelet_initialization.py @@ -0,0 +1,144 @@ +#!/usr/bin/env python + +## \file run.py +# \brief Flame simulation initialization test case. +# \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 . + +import pysu2 +from mpi4py import MPI + +# Generatl config options +config_options = """ +FILENAMES_INTERPOLATOR= (MLP_PIML_1.mlp,MLP_PIML_2.mlp,MLP_PIML_3.mlp,MLP_PIML_4.mlp,MLP_PIML_5.mlp,MLP_PIML_6.mlp,MLP_PIML_7.mlp,MLP_PIML_8.mlp,MLP_PIML_9.mlp,MLP_PIML_10.mlp,MLP_PIML_11.mlp) + +INC_VELOCITY_INIT= (1e-5, 0.0, 0.0) +MARKER_SYM=(x_minus, x_plus, y_plus, y_minus) + +SOLVER = INC_NAVIER_STOKES +INC_DENSITY_MODEL= VARIABLE +INC_ENERGY_EQUATION = YES +INC_TEMPERATURE_INIT= 300 +INC_NONDIM= DIMENSIONAL + + +FLUID_MODEL= FLUID_FLAMELET +INTERPOLATION_METHOD= MLP + +CONTROLLING_VARIABLE_NAMES=(ProgressVariable,EnthalpyTot,MixtureFraction) +CONTROLLING_VARIABLE_SOURCE_NAMES=(ProdRateTot_PV,NULL,NULL) + +KIND_SCALAR_MODEL= FLAMELET +DIFFUSIVITY_MODEL= FLAMELET +VISCOSITY_MODEL= FLAMELET +CONDUCTIVITY_MODEL= FLAMELET + +PREFERENTIAL_DIFFUSION=YES + +FLAME_INIT_METHOD= __FLAME_INIT_METHOD__ +SPARK_INIT=(5e-4, 2.5e-4, 0.0, 1e-4, 1, 50) +SPARK_REACTION_RATES=(1.5e4, 0,0) +FLAME_INIT= (4e-4, 0.00, 0.00, 1.0, 0.0, 0.0, 5e-5, 1.0) +FLAME_INIT_IGNITION= __IGNITION_TEMP__ + +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR + +MUSCL_SPECIES= YES +TIME_DISCRE_SPECIES= EULER_IMPLICIT +SPECIES_CLIPPING= NO +SPECIES_INIT=(-2.251135e-01,+2.262913e+03,+1.446752e-02) + +CFL_REDUCTION_SPECIES= 1.0 + +INC_OUTLET_TYPE= PRESSURE_OUTLET + +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 50 +CFL_ADAPT= NO +ITER= __N_ITER__ + +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 1E-03 +LINEAR_SOLVER_ITER= 15 + +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +TIME_DISCRE_FLOW= EULER_IMPLICIT + +SCREEN_OUTPUT = CUSTOM +HISTORY_OUTPUT = CUSTOM +VOLUME_OUTPUT = SOLUTION + +MESH_FORMAT= RECTANGLE +MESH_BOX_SIZE= ( 200, 50, 0 ) +MESH_BOX_LENGTH= ( 1e-3, 5e-4, 0 ) + +OUTPUT_FILES=(RESTART) +CUSTOM_OUTPUTS= 'Tprobe : Probe{TEMPERATURE}[5e-4, 2.5e-4];' + +""" + +comm = MPI.COMM_WORLD + +def main(): + + comm = MPI.COMM_WORLD + + T_ignition = 1600 + if comm.Get_rank() == 0: + + # Prepare configuration files for flame front and spark ignition methods. + config_flame_front = config_options.replace("__FLAME_INIT_METHOD__", "FLAME_FRONT") + config_flame_front = config_flame_front.replace("__N_ITER__", "1") + config_flame_front = config_flame_front.replace("__IGNITION_TEMP__", "%f" % T_ignition) + + config_spark = config_options.replace("__FLAME_INIT_METHOD__", "SPARK") + config_spark = config_spark.replace("__N_ITER__", "20") + config_spark = config_spark.replace("__IGNITION_TEMP__", "%f" % T_ignition) + + with open("config_flame_front.cfg", "w") as f: + f.write(config_flame_front) + with open("config_spark.cfg", "w") as f: + f.write(config_spark) + comm.Barrier() + + # Initialize simulation with flame front ignition + driver_flame_front = pysu2.CSinglezoneDriver("config_flame_front.cfg", 1, comm) + driver_flame_front.StartSolver() + + # Retrieve temperature from the burnt zone of the initial flame front. + T_flame_front = driver_flame_front.GetOutputValue("Tprobe") + driver_flame_front.Finalize() + + # Initialize simulation with spark ignition. + driver_spark = pysu2.CSinglezoneDriver("config_spark.cfg", 1, comm) + driver_spark.StartSolver() + + # Retrieve temperature from center of the spark. + T_spark = driver_spark.GetOutputValue("Tprobe") + driver_spark.Finalize() + if comm.Get_rank() == 0: + print("| %i | %.2f | %.2f|" % (1, min(T_ignition, T_flame_front), min(T_ignition, T_spark))) + +if __name__ == '__main__': + main() diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 1c25e4942d22..9b5ebf1344a3 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -75,9 +75,17 @@ def main(): cfd_flamelet_h2.cfg_file = "laminar_premixed_h2_flame_cfd.cfg" cfd_flamelet_h2.test_iter = 5 cfd_flamelet_h2.test_vals = [-8.036958, -8.372668, -1.842800, -9.388446] - cfd_flamelet_h2.new_output = True test_list.append(cfd_flamelet_h2) + # Flame ignition methods + flame_init_methods = TestCase("flame_init_methods") + flame_init_methods.cfg_dir = "flamelet/08_flame_init_method" + flame_init_methods.cfg_file ="flamelet_initialization.py" + flame_init_methods.command = TestCase.Command("mpirun -np 2", "python", "flamelet_initialization.py") + flame_init_methods.test_vals=[1600, 1600] + flame_init_methods.new_output = True + test_list.append(flame_init_methods) + ######################### ## NEMO solver ### ######################### diff --git a/config_template.cfg b/config_template.cfg index ab8175bb6802..f5fbdeb96dc2 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -960,6 +960,16 @@ PREFERENTIAL_DIFFUSION= NO % (x8) = Thickness of the 'burnt' zone, after this length, the conditions will be 'unburnt' again. FLAME_INIT= (0.004, 0.0, 0.0, 1.0, 0.0, 0.0, 0.2e-3, 1.0) +% Flame ignition temperature in Kelvin. +% Specify the temperature on the 'burnt' side of the flame when using FLAME_FRONT initialization +% or the maximum temperature within the spark region when using SPARK initialization. +% The initial value of the progress variable on the 'burnt' side of the flame front is calculated through +% inverse regression of the ignition temperature or when reaching the upper limit of the progress variable +% in the flamelet data space. +% The SPARK_REACTION_RATES are applied within the spark region when the local temperature is below the +% flame ignition temperature. +FLAME_INIT_IGNITION= 5000.0 + % SPARK initialization % the solution is ignited through application of a set of source terms within a specified region for a set number % of solver iterations. This artificial spark can be applied after a certain number of iterations has passed, allowing for diff --git a/meson_scripts/init.py b/meson_scripts/init.py index 858456015f13..9e8d18dcde75 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 = "02f2cb9dde791074858e11ac091f7c4df9c6af65" + sha_version_mlpcpp = "e23facf388902f262fbe7ba3bcc84d36c85350b9" 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 02f2cb9dde79..e23facf38890 160000 --- a/subprojects/MLPCpp +++ b/subprojects/MLPCpp @@ -1 +1 @@ -Subproject commit 02f2cb9dde791074858e11ac091f7c4df9c6af65 +Subproject commit e23facf388902f262fbe7ba3bcc84d36c85350b9