From c80aa76995c46d6b3bd2498693fb39986e453594 Mon Sep 17 00:00:00 2001 From: capriatim Date: Fri, 6 Jun 2025 11:51:20 +0200 Subject: [PATCH 1/5] expanded B' to specify condesed phase --- interface/fortran/cwrapper.cpp | 11 ++++++-- src/apps/bprime.cpp | 41 +++++++++++++++++++++------- src/thermo/Thermodynamics.cpp | 49 ++++++++++++++++++++++------------ src/thermo/Thermodynamics.h | 5 ++-- 4 files changed, 75 insertions(+), 31 deletions(-) diff --git a/interface/fortran/cwrapper.cpp b/interface/fortran/cwrapper.cpp index 777a78ab..dac4100b 100644 --- a/interface/fortran/cwrapper.cpp +++ b/interface/fortran/cwrapper.cpp @@ -452,7 +452,11 @@ void NAME_MANGLE(surface_mass_balance) const double* const P, const double* const Bg, double* const Bc, double* const hw, double *const p_Xs) { - p_mix->surfaceMassBalance(p_Yke, p_Ykg, *T, *P, *Bg, *Bc, *hw, p_Xs); + const int ne = p_mix->nElements(); + std::vector Ykc (ne,0); + int ic = p_mix->elementIndex("C"); + Ykc[ic] = 1.0; + p_mix->surfaceMassBalance(p_Yke, p_Ykg, Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); } //============================================================================== @@ -474,11 +478,14 @@ void NAME_MANGLE(gasmixture_surface_mass_balance) std::vector Yke (ne,0); std::vector Ykg (ne,0); + std::vector Ykc (ne,0); + int ic = p_mix->elementIndex("C"); + Ykc[ic] = 1.0; p_mix->getComposition(char_to_string(edge, edge_length), Yke.data(), Composition::MASS); p_mix->getComposition(char_to_string(pyro, pyro_length), Ykg.data(), Composition::MASS); - p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); + p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); } //============================================================================== diff --git a/src/apps/bprime.cpp b/src/apps/bprime.cpp index 3de7ec24..6a246b0d 100644 --- a/src/apps/bprime.cpp +++ b/src/apps/bprime.cpp @@ -46,17 +46,16 @@ using namespace Mutation::Utilities; * * bprime \f$-T\f$ \f$T_1\f$:\f$\Delta T\f$:\f$T_2\f$ \f$-p\f * \f$p\f$ \f$-b\f \f$B'_g\f$ \f$-m\f mixture \f$-bl\f BL \f$-py \f Pyrolysis - * + * \f$-cp \f CondencedPhase * This program generates a so-called "B-prime" table for a given temperature * range and stepsize in K, a fixed pressure in Pa, a value of \f$B'_g\f$ * (pyrolysis mass flux, nondimensionalized by the boundary layer edge mass - * flux), a given mixture name. Currently, this program assumes the char - * composition to be solid graphite (which must be included in the mixture - * file). The `BL` and `Pyrolysis` arguments are the names of the [elemental - * compositions](@ref compositions) in the mixture file which represent the - * boundary layer edge and pyrolysis gases respectively. The produced table - * provides values of \f$B'_c\f$, the wall enthalpy in MJ/kg, and the species - * mole fractions at the wall versus temperature. + * flux), a given mixture name. + * The `BL`, `Pyrolysis`, and `CondencedPhase` arguments are the names of the + * [elemental compositions](@ref compositions) in the mixture file which represent + * the boundary layer edge and pyrolysis gases, and the condenced phase. + * The produced table provides values of \f$B'_c\f$, the wall enthalpy in MJ/kg, + * and the species mole fractions at the wall versus temperature. */ // Simply stores the command line options typedef struct Options { @@ -70,8 +69,10 @@ typedef struct Options { std::string mixture; std::string boundary_layer_comp; std::string pyrolysis_composition; + std::string condenced_phase_composition; bool pyrolysis_exist = false; + bool condenced_phase_exist = false; } Options; // Checks if an option is present @@ -122,12 +123,15 @@ void printHelpMessage(const char* const name) { cout << tab << "-py pyrolysis composition name (default = null)" << endl; + cout << tab + << "-cp condenced phase composition name (default = carbon)" + << endl; cout << endl; cout << "Example:" << endl; cout << tab << name - << " -T 300:100:5000 -P 101325 -b 10 -m carbonPhenol -bl BLedge -py Gas" + << " -T 300:100:5000 -P 101325 -b 10 -m carbonPhenol -bl BLedge -py Gas -cp CondencedPhase" << endl; cout << endl; cout << "Mixture file:" << endl; @@ -140,6 +144,9 @@ void printHelpMessage(const char* const name) { cout << tab << "Gas - corresponds to the pyrolysis elemental gas composition" << endl; + cout << tab + << "CondencedPhase - corresponds to the condenced phase composition" + << endl; cout << endl; exit(0); @@ -246,6 +253,12 @@ Options parseOptions(int argc, char** argv) { opts.pyrolysis_exist = true; } + if (optionExists(argc, argv, "-cp")) { + opts.condenced_phase_composition = getOption(argc, argv, "-cp"); + opts.condenced_phase_exist = true; + } + + return opts; } @@ -263,6 +276,7 @@ int main(int argc, char* argv[]) { std::vector Yke(ne, 0); std::vector Ykg(ne, 0); + std::vector Ycp(ne, 0); std::vector Xw(ns, 0); // Run conditions @@ -279,6 +293,13 @@ int main(int argc, char* argv[]) { mix.getComposition(opts.pyrolysis_composition, Ykg.data(), Composition::MASS); + if (opts.condenced_phase_exist) + mix.getComposition(opts.condenced_phase_composition, Ycp.data(), Composition::MASS); + else { + int ic = mix.elementIndex("C"); + Ycp[ic] = 1.0; + } + cout << setw(10) << "\"Tw[K]\"" << setw(15) << "\"B'c\"" << setw(15) << "\"hw[MJ/kg]\""; for (int i = 0; i < ns; ++i) @@ -286,7 +307,7 @@ int main(int argc, char* argv[]) { cout << endl; for (double T = T1; T < T2 + 1.0e-6; T += dt) { - mix.surfaceMassBalance(Yke.data(), Ykg.data(), T, P, Bg, Bc, hw, + mix.surfaceMassBalance(Yke.data(), Ykg.data(), Ycp.data(), T, P, Bg, Bc, hw, Xw.data()); cout << setw(10) << T << setw(15) << Bc << setw(15) << hw / 1.0e6; for (int i = 0; i < ns; ++i) cout << setw(25) << Xw[i]; diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index ddad37ac..78fe963b 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -1090,8 +1090,9 @@ void Thermodynamics::elementFractions( //============================================================================== void Thermodynamics::surfaceMassBalance( - const double *const p_Yke, const double *const p_Ykg, const double T, - const double P, const double Bg, double &Bc, double &hw, double *const p_Xs) + const double *const p_Yke, const double *const p_Ykg, const double *const p_Ycp, + const double T, const double P, const double Bg, double &Bc, double &hw, + double *const p_Xs) { const int ne = nElements(); const int ng = nGas(); @@ -1107,12 +1108,16 @@ void Thermodynamics::surfaceMassBalance( sum += p_Xw[i]; } - // Use "large" amount of carbon to simulate infinite char - int ic = elementIndex("C"); - //double carbon = std::min(1000.0, std::max(100.0,1000.0*Bg)); - double carbon = std::max(100.0*Bg, 200.0); - p_Xw[ic] += carbon; - sum += carbon; + // Use "large" amount of condences phase to simulate infinite surface + double LargeNumber = 100.0; + std::vector condencedPhaseElements; + double tol = 1.0e-16; + for (int i = 0; i < ne; ++i) { + p_Xw[i] += LargeNumber*p_Ycp[i]; + sum += LargeNumber*p_Ycp[i]; + if (abs(p_Ycp[i]) > tol) + condencedPhaseElements.push_back(i); + } for (int i = 0; i < ne; ++i) p_Xw[i] /= sum; @@ -1123,22 +1128,32 @@ void Thermodynamics::surfaceMassBalance( // Compute the gas mass fractions at the wall double mwg = 0.0; - double ywc = 0.0; + double p_Yw[ne]; for (int j = 0; j < ng; ++j) { mwg += speciesMw(j) * p_X[j]; - ywc += elementMatrix()(j,ic) * p_X[j]; - //for (int i = 0; i < ne; ++i) - // p_Yw[i] += elementMatrix()(j,i) * p_X[j]; + for (int i = 0; i < ne; ++i) + p_Yw[i] += elementMatrix()(j,i) * p_X[j]; } - //for (int i = 0; i < ne; ++i) - // p_Yw[i] *= atomicMass(i) / mwg; - ywc *= atomicMass(ic) / mwg; + for (int i = 0; i < ne; ++i) + p_Yw[i] *= atomicMass(i) / mwg; + + double sum_Ye = 0.0; + double sum_Yg = 0.0; + double sum_Yw = 0.0; + double sum_YCp = 1.0; //Assumed equal to one + int ncp = condencedPhaseElements.size(); + for (int i=0; i < ncp; ++i ) { + int idx_cp = condencedPhaseElements[i]; + sum_Ye += p_Yke[idx_cp]; + sum_Yg += p_Ykg[idx_cp]; + sum_Yw += p_Yw[idx_cp]; + } // Compute char mass blowing rate - Bc = (p_Yke[ic] + Bg*p_Ykg[ic] - ywc*(1.0 + Bg)) / (ywc - 1.0); - Bc = std::max(Bc, 0.0); + Bc = (Bg*(sum_Yg - sum_Yw) + sum_Ye - sum_Yw)/(sum_Yw - sum_YCp); + Bc = std::max(Bc, 1.0e-8); // Compute the gas enthalpy speciesHOverRT(T, p_h); diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index 9ba1b888..e05e723f 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -873,6 +873,7 @@ class Thermodynamics //: public StateModelUpdateHandler * * @param p_Yke Element mass fractions of boundary layer edge. * @param p_Ykg Element mass fractions of the pyrolysis gas. + * @param p_Ycp Element mass fractions of the condenced phase. * @param T Temperature at the surface in K. * @param P Pressure at the surface in Pa. * @param Bg Non-dimensional pyrolysis gas mass blowing rate. @@ -882,8 +883,8 @@ class Thermodynamics //: public StateModelUpdateHandler * the surface. */ void surfaceMassBalance( - const double *const p_Yke, const double *const p_Ykg, const double T, - const double P, const double Bg, double &Bc, double &hw, + const double *const p_Yke, const double *const p_Ykg, const double *const p_Ycp, + const double T, const double P, const double Bg, double &Bc, double &hw, double *const p_Xs = NULL); private: From 272a48fe46d8c81c0d240d7f042572d2c4a44843 Mon Sep 17 00:00:00 2001 From: capriatim Date: Tue, 9 Jun 2026 11:42:08 +0200 Subject: [PATCH 2/5] Addressed Feedbacks Pull Request #279 --- interface/fortran/cwrapper.cpp | 29 +++++++++++++++++++++++++++++ interface/fortran/cwrapper.h | 22 ++++++++++++++++++++-- src/apps/bprime.cpp | 30 +++++++++++++++--------------- src/thermo/Thermodynamics.cpp | 18 +++++++++--------- src/thermo/Thermodynamics.h | 4 ++-- 5 files changed, 75 insertions(+), 28 deletions(-) diff --git a/interface/fortran/cwrapper.cpp b/interface/fortran/cwrapper.cpp index dac4100b..8674449c 100644 --- a/interface/fortran/cwrapper.cpp +++ b/interface/fortran/cwrapper.cpp @@ -459,6 +459,15 @@ void NAME_MANGLE(surface_mass_balance) p_mix->surfaceMassBalance(p_Yke, p_Ykg, Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); } +//============================================================================== +void NAME_MANGLE(surface_mass_balance_general) + (const double *const p_Yke, const double *const p_Ykg, const double *const p_Ykc, + const double* const T, const double* const P, const double* const Bg, + double* const Bc, double* const hw, double *const p_Xs) +{ + p_mix->surfaceMassBalance(p_Yke, p_Ykg, p_Ykc, *T, *P, *Bg, *Bc, *hw, p_Xs); +} + //============================================================================== void NAME_MANGLE(get_composition) (F_STRING mixture, double* const p_Yk, F_STRLEN mixture_length) @@ -488,6 +497,26 @@ void NAME_MANGLE(gasmixture_surface_mass_balance) p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); } +//============================================================================== +void NAME_MANGLE(gasmixture_surface_mass_balance_general) + (F_STRING edge, F_STRING pyro, F_STRING surf, + const double* const T, const double* const P, const double* const Bg, + double* const Bc, double* const hw, double *const p_Xs, + F_STRLEN edge_length, F_STRLEN pyro_length, F_STRLEN surf_length) +{ + const int ne = p_mix->nElements(); + + std::vector Yke (ne,0); + std::vector Ykg (ne,0); + std::vector Ykc (ne,0); + + p_mix->getComposition(char_to_string(edge, edge_length), Yke.data(), Composition::MASS); + p_mix->getComposition(char_to_string(pyro, pyro_length), Ykg.data(), Composition::MASS); + p_mix->getComposition(char_to_string(surf, surf_length), Ykc.data(), Composition::MASS); + + p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); +} + //============================================================================== void NAME_MANGLE(source_energy_transfer) (double *const p_source_transfer) diff --git a/interface/fortran/cwrapper.h b/interface/fortran/cwrapper.h index 0358caa9..287e894c 100644 --- a/interface/fortran/cwrapper.h +++ b/interface/fortran/cwrapper.h @@ -444,13 +444,22 @@ void NAME_MANGLE(diffusion_matrix)(double* const p_Dij); double NAME_MANGLE(sigma)(); /** - * Solves the surface mass balance at an ablating surface. + * Solves the surface mass balance at a carbon ablating surface. */ void NAME_MANGLE(surface_mass_balance) (const double *const p_Yke, const double *const p_Ykg, const double* const T, const double* const P, const double* const Bg, double* const Bc, double* const hw, double *const p_Xs); +/** + * Solves the surface mass balance at a general ablating surface. + */ +void NAME_MANGLE(surface_mass_balance_general) + (const double *const p_Yke, const double *const p_Ykg, const double *const p_Ykc, + const double* const T, const double* const P, const double* const Bg, + double* const Bc, double* const hw, double *const p_Xs); + + /** * Gets the composition for components of the mixture */ @@ -458,7 +467,7 @@ void NAME_MANGLE(get_composition) (F_STRING mixture, double* const p_Yk, F_STRLEN mixture_length); /** - * Solves the surface mass balance at an ablating surface provided gas mixture names. + * Solves the surface mass balance at a carbon ablating surface provided gas mixture names. */ void NAME_MANGLE(gasmixture_surface_mass_balance) (F_STRING edge, F_STRING pyro, @@ -466,6 +475,15 @@ void NAME_MANGLE(gasmixture_surface_mass_balance) double* const Bc, double* const hw, double *const p_Xs, F_STRLEN edge_length, F_STRLEN pyro_length); +/** + * Solves the surface mass balance at a general ablating surface provided gas mixture names. + */ +void NAME_MANGLE(gasmixture_surface_mass_balance_general) + (F_STRING edge, F_STRING pyro, F_STRING surf, + const double* const T, const double* const P, const double* const Bg, + double* const Bc, double* const hw, double *const p_Xs, + F_STRLEN edge_length, F_STRLEN pyro_length, F_STRLEN surf_length); + /** * Returns the pointer to the energy transfer between the internal * energy modes. diff --git a/src/apps/bprime.cpp b/src/apps/bprime.cpp index 6a246b0d..141c2aba 100644 --- a/src/apps/bprime.cpp +++ b/src/apps/bprime.cpp @@ -46,14 +46,14 @@ using namespace Mutation::Utilities; * * bprime \f$-T\f$ \f$T_1\f$:\f$\Delta T\f$:\f$T_2\f$ \f$-p\f * \f$p\f$ \f$-b\f \f$B'_g\f$ \f$-m\f mixture \f$-bl\f BL \f$-py \f Pyrolysis - * \f$-cp \f CondencedPhase + * \f$-cp \f CondensedPhase * This program generates a so-called "B-prime" table for a given temperature * range and stepsize in K, a fixed pressure in Pa, a value of \f$B'_g\f$ * (pyrolysis mass flux, nondimensionalized by the boundary layer edge mass * flux), a given mixture name. - * The `BL`, `Pyrolysis`, and `CondencedPhase` arguments are the names of the + * The `BL`, `Pyrolysis`, and `CondensedPhase` arguments are the names of the * [elemental compositions](@ref compositions) in the mixture file which represent - * the boundary layer edge and pyrolysis gases, and the condenced phase. + * the boundary layer edge and pyrolysis gases, and the condensed phase. * The produced table provides values of \f$B'_c\f$, the wall enthalpy in MJ/kg, * and the species mole fractions at the wall versus temperature. */ @@ -69,10 +69,10 @@ typedef struct Options { std::string mixture; std::string boundary_layer_comp; std::string pyrolysis_composition; - std::string condenced_phase_composition; + std::string condensed_phase_composition; bool pyrolysis_exist = false; - bool condenced_phase_exist = false; + bool condensed_phase_exist = false; } Options; // Checks if an option is present @@ -124,14 +124,14 @@ void printHelpMessage(const char* const name) { << "-py pyrolysis composition name (default = null)" << endl; cout << tab - << "-cp condenced phase composition name (default = carbon)" + << "-cp condensed phase composition name (default = carbon)" << endl; cout << endl; cout << "Example:" << endl; cout << tab << name - << " -T 300:100:5000 -P 101325 -b 10 -m carbonPhenol -bl BLedge -py Gas -cp CondencedPhase" + << " -T 300:100:5000 -P 101325 -b 10 -m carbonPhenol -bl BLedge -py Gas -cp CondensedPhase" << endl; cout << endl; cout << "Mixture file:" << endl; @@ -145,7 +145,7 @@ void printHelpMessage(const char* const name) { << "Gas - corresponds to the pyrolysis elemental gas composition" << endl; cout << tab - << "CondencedPhase - corresponds to the condenced phase composition" + << "CondensedPhase - corresponds to the condensed phase composition" << endl; cout << endl; @@ -254,8 +254,8 @@ Options parseOptions(int argc, char** argv) { } if (optionExists(argc, argv, "-cp")) { - opts.condenced_phase_composition = getOption(argc, argv, "-cp"); - opts.condenced_phase_exist = true; + opts.condensed_phase_composition = getOption(argc, argv, "-cp"); + opts.condensed_phase_exist = true; } @@ -276,7 +276,7 @@ int main(int argc, char* argv[]) { std::vector Yke(ne, 0); std::vector Ykg(ne, 0); - std::vector Ycp(ne, 0); + std::vector Ykc(ne, 0); std::vector Xw(ns, 0); // Run conditions @@ -293,11 +293,11 @@ int main(int argc, char* argv[]) { mix.getComposition(opts.pyrolysis_composition, Ykg.data(), Composition::MASS); - if (opts.condenced_phase_exist) - mix.getComposition(opts.condenced_phase_composition, Ycp.data(), Composition::MASS); + if (opts.condensed_phase_exist) + mix.getComposition(opts.condensed_phase_composition, Ykc.data(), Composition::MASS); else { int ic = mix.elementIndex("C"); - Ycp[ic] = 1.0; + Ykc[ic] = 1.0; } cout << setw(10) << "\"Tw[K]\"" << setw(15) << "\"B'c\"" << setw(15) @@ -307,7 +307,7 @@ int main(int argc, char* argv[]) { cout << endl; for (double T = T1; T < T2 + 1.0e-6; T += dt) { - mix.surfaceMassBalance(Yke.data(), Ykg.data(), Ycp.data(), T, P, Bg, Bc, hw, + mix.surfaceMassBalance(Yke.data(), Ykg.data(), Ykc.data(), T, P, Bg, Bc, hw, Xw.data()); cout << setw(10) << T << setw(15) << Bc << setw(15) << hw / 1.0e6; for (int i = 0; i < ns; ++i) cout << setw(25) << Xw[i]; diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index 78fe963b..8006b64d 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -1090,7 +1090,7 @@ void Thermodynamics::elementFractions( //============================================================================== void Thermodynamics::surfaceMassBalance( - const double *const p_Yke, const double *const p_Ykg, const double *const p_Ycp, + const double *const p_Yke, const double *const p_Ykg, const double *const p_Ykc, const double T, const double P, const double Bg, double &Bc, double &hw, double *const p_Xs) { @@ -1110,13 +1110,13 @@ void Thermodynamics::surfaceMassBalance( // Use "large" amount of condences phase to simulate infinite surface double LargeNumber = 100.0; - std::vector condencedPhaseElements; + std::vector condensedPhaseElements; double tol = 1.0e-16; for (int i = 0; i < ne; ++i) { - p_Xw[i] += LargeNumber*p_Ycp[i]; - sum += LargeNumber*p_Ycp[i]; - if (abs(p_Ycp[i]) > tol) - condencedPhaseElements.push_back(i); + p_Xw[i] += LargeNumber*p_Ykc[i]; + sum += LargeNumber*p_Ykc[i]; + if (abs(p_Ykc[i]) > tol) + condensedPhaseElements.push_back(i); } for (int i = 0; i < ne; ++i) @@ -1143,9 +1143,9 @@ void Thermodynamics::surfaceMassBalance( double sum_Yg = 0.0; double sum_Yw = 0.0; double sum_YCp = 1.0; //Assumed equal to one - int ncp = condencedPhaseElements.size(); + int ncp = condensedPhaseElements.size(); for (int i=0; i < ncp; ++i ) { - int idx_cp = condencedPhaseElements[i]; + int idx_cp = condensedPhaseElements[i]; sum_Ye += p_Yke[idx_cp]; sum_Yg += p_Ykg[idx_cp]; sum_Yw += p_Yw[idx_cp]; @@ -1153,7 +1153,7 @@ void Thermodynamics::surfaceMassBalance( // Compute char mass blowing rate Bc = (Bg*(sum_Yg - sum_Yw) + sum_Ye - sum_Yw)/(sum_Yw - sum_YCp); - Bc = std::max(Bc, 1.0e-8); + Bc = std::max(Bc, 1.0e-16); // Compute the gas enthalpy speciesHOverRT(T, p_h); diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index e05e723f..b457c6d9 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -873,7 +873,7 @@ class Thermodynamics //: public StateModelUpdateHandler * * @param p_Yke Element mass fractions of boundary layer edge. * @param p_Ykg Element mass fractions of the pyrolysis gas. - * @param p_Ycp Element mass fractions of the condenced phase. + * @param p_Ykc Element mass fractions of the condensed phase. * @param T Temperature at the surface in K. * @param P Pressure at the surface in Pa. * @param Bg Non-dimensional pyrolysis gas mass blowing rate. @@ -883,7 +883,7 @@ class Thermodynamics //: public StateModelUpdateHandler * the surface. */ void surfaceMassBalance( - const double *const p_Yke, const double *const p_Ykg, const double *const p_Ycp, + const double *const p_Yke, const double *const p_Ykg, const double *const p_Ykc, const double T, const double P, const double Bg, double &Bc, double &hw, double *const p_Xs = NULL); From 5425e148be479e832ee95e37a6f55370cf2f6684 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 10 Jun 2026 10:43:28 +0200 Subject: [PATCH 3/5] Addressed Feedbacks Pull Request mutationpp#279 ensuring backward-compatibility --- CMakeLists.txt | 2 +- interface/fortran/cwrapper.cpp | 15 ++++----------- src/apps/bprime.cpp | 4 ++-- src/thermo/Thermodynamics.cpp | 26 ++++++++++++++++++++------ src/thermo/Thermodynamics.h | 7 ++++--- 5 files changed, 31 insertions(+), 23 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ff9ec2b..6c0242f7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,7 +32,7 @@ cmake_policy(SET CMP0048 NEW) set(CMAKE_CXX_STANDARD 14) project(mutation++ - VERSION 1.1.3 + VERSION 1.2.0 LANGUAGES CXX ) diff --git a/interface/fortran/cwrapper.cpp b/interface/fortran/cwrapper.cpp index 8674449c..4c19080e 100644 --- a/interface/fortran/cwrapper.cpp +++ b/interface/fortran/cwrapper.cpp @@ -452,11 +452,7 @@ void NAME_MANGLE(surface_mass_balance) const double* const P, const double* const Bg, double* const Bc, double* const hw, double *const p_Xs) { - const int ne = p_mix->nElements(); - std::vector Ykc (ne,0); - int ic = p_mix->elementIndex("C"); - Ykc[ic] = 1.0; - p_mix->surfaceMassBalance(p_Yke, p_Ykg, Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); + p_mix->surfaceMassBalance(p_Yke, p_Ykg, *T, *P, *Bg, *Bc, *hw, p_Xs); } //============================================================================== @@ -465,7 +461,7 @@ void NAME_MANGLE(surface_mass_balance_general) const double* const T, const double* const P, const double* const Bg, double* const Bc, double* const hw, double *const p_Xs) { - p_mix->surfaceMassBalance(p_Yke, p_Ykg, p_Ykc, *T, *P, *Bg, *Bc, *hw, p_Xs); + p_mix->surfaceMassBalance(p_Yke, p_Ykg, *T, *P, *Bg, *Bc, *hw, p_Xs, p_Ykc); } //============================================================================== @@ -487,14 +483,11 @@ void NAME_MANGLE(gasmixture_surface_mass_balance) std::vector Yke (ne,0); std::vector Ykg (ne,0); - std::vector Ykc (ne,0); - int ic = p_mix->elementIndex("C"); - Ykc[ic] = 1.0; p_mix->getComposition(char_to_string(edge, edge_length), Yke.data(), Composition::MASS); p_mix->getComposition(char_to_string(pyro, pyro_length), Ykg.data(), Composition::MASS); - p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); + p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); } //============================================================================== @@ -514,7 +507,7 @@ void NAME_MANGLE(gasmixture_surface_mass_balance_general) p_mix->getComposition(char_to_string(pyro, pyro_length), Ykg.data(), Composition::MASS); p_mix->getComposition(char_to_string(surf, surf_length), Ykc.data(), Composition::MASS); - p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs); + p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), *T, *P, *Bg, *Bc, *hw, p_Xs, Ykc.data()); } //============================================================================== diff --git a/src/apps/bprime.cpp b/src/apps/bprime.cpp index 141c2aba..73832796 100644 --- a/src/apps/bprime.cpp +++ b/src/apps/bprime.cpp @@ -307,8 +307,8 @@ int main(int argc, char* argv[]) { cout << endl; for (double T = T1; T < T2 + 1.0e-6; T += dt) { - mix.surfaceMassBalance(Yke.data(), Ykg.data(), Ykc.data(), T, P, Bg, Bc, hw, - Xw.data()); + mix.surfaceMassBalance(Yke.data(), Ykg.data(), T, P, Bg, Bc, hw, + Xw.data(), Ykc.data()); cout << setw(10) << T << setw(15) << Bc << setw(15) << hw / 1.0e6; for (int i = 0; i < ns; ++i) cout << setw(25) << Xw[i]; cout << endl; diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index 8006b64d..3d3b5de7 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -1090,9 +1090,9 @@ void Thermodynamics::elementFractions( //============================================================================== void Thermodynamics::surfaceMassBalance( - const double *const p_Yke, const double *const p_Ykg, const double *const p_Ykc, + const double *const p_Yke, const double *const p_Ykg, const double T, const double P, const double Bg, double &Bc, double &hw, - double *const p_Xs) + double *const p_Xs, const double *const p_Ykc) { const int ne = nElements(); const int ng = nGas(); @@ -1101,6 +1101,17 @@ void Thermodynamics::surfaceMassBalance( double* p_X = (p_Xs != NULL ? p_Xs : mp_work1); double* p_h = mp_work2; + const double *p_Ykc_to_use = p_Ykc; + std::vector default_Ykc(ne); + if (p_Ykc == NULL) { + for (int i = 0; i < ne; ++i) { + default_Ykc[i] = 0.0; + } + int ic = elementIndex("C"); + default_Ykc[ic] = 1.0; + p_Ykc_to_use = default_Ykc.data(); + } + // Initialize the wall element fractions to be the pyrolysis gas fractions double sum = 0.0; for (int i = 0; i < ne; ++i) { @@ -1113,9 +1124,9 @@ void Thermodynamics::surfaceMassBalance( std::vector condensedPhaseElements; double tol = 1.0e-16; for (int i = 0; i < ne; ++i) { - p_Xw[i] += LargeNumber*p_Ykc[i]; - sum += LargeNumber*p_Ykc[i]; - if (abs(p_Ykc[i]) > tol) + p_Xw[i] += LargeNumber*p_Ykc_to_use[i]; + sum += LargeNumber*p_Ykc_to_use[i]; + if (abs(p_Ykc_to_use[i]) > tol) condensedPhaseElements.push_back(i); } @@ -1129,6 +1140,9 @@ void Thermodynamics::surfaceMassBalance( // Compute the gas mass fractions at the wall double mwg = 0.0; double p_Yw[ne]; + for (int i = 0; i < ne; ++i) { + p_Yw[i] = 0.0; + } for (int j = 0; j < ng; ++j) { mwg += speciesMw(j) * p_X[j]; @@ -1153,7 +1167,7 @@ void Thermodynamics::surfaceMassBalance( // Compute char mass blowing rate Bc = (Bg*(sum_Yg - sum_Yw) + sum_Ye - sum_Yw)/(sum_Yw - sum_YCp); - Bc = std::max(Bc, 1.0e-16); + Bc = std::max(Bc, 1.0e-15); // Compute the gas enthalpy speciesHOverRT(T, p_h); diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index b457c6d9..92f7df64 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -873,7 +873,6 @@ class Thermodynamics //: public StateModelUpdateHandler * * @param p_Yke Element mass fractions of boundary layer edge. * @param p_Ykg Element mass fractions of the pyrolysis gas. - * @param p_Ykc Element mass fractions of the condensed phase. * @param T Temperature at the surface in K. * @param P Pressure at the surface in Pa. * @param Bg Non-dimensional pyrolysis gas mass blowing rate. @@ -881,11 +880,13 @@ class Thermodynamics //: public StateModelUpdateHandler * @param hw On return, enthalpy of the mixed gas at the surface in J/kg. * @param p_Xs (optional) On return, vector of species mole fractions at * the surface. + * @param p_Ykc (optional) Element mass fractions of the condensed phase. + * If not value is passed, a carbon surface is assumed. */ void surfaceMassBalance( - const double *const p_Yke, const double *const p_Ykg, const double *const p_Ykc, + const double *const p_Yke, const double *const p_Ykg, const double T, const double P, const double Bg, double &Bc, double &hw, - double *const p_Xs = NULL); + double *const p_Xs = NULL, const double *const p_Ykc = NULL); private: From a630a5a14948d2757265687ff9084015f9fa2ed0 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 10 Jun 2026 10:44:37 +0200 Subject: [PATCH 4/5] Solving issue #299 --- examples/fortran/wrapper_test/wrapper_test.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fortran/wrapper_test/wrapper_test.f90 b/examples/fortran/wrapper_test/wrapper_test.f90 index 3a7361f8..458e66cf 100644 --- a/examples/fortran/wrapper_test/wrapper_test.f90 +++ b/examples/fortran/wrapper_test/wrapper_test.f90 @@ -74,7 +74,7 @@ program main do i = 1,295 T = dble(i-1)*50.0 + 300.0 - call mpp_set_state(T, P, var) + call mpp_set_state(P, T, var) call mpp_x(species_x) write(*,'(E12.4)',advance='no') T From a60b6f10af536b5c4dbc59985e2a9b175902ed62 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 10 Jun 2026 16:51:05 +0200 Subject: [PATCH 5/5] Addressed Feedbacks Pull Request mutationpp#279 pt3 --- src/thermo/Thermodynamics.cpp | 10 +++++----- src/thermo/Thermodynamics.h | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index 3d3b5de7..7febb2c8 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -1098,12 +1098,12 @@ void Thermodynamics::surfaceMassBalance( const int ng = nGas(); double p_Xw [ne]; - double* p_X = (p_Xs != NULL ? p_Xs : mp_work1); + double* p_X = (p_Xs != nullptr ? p_Xs : mp_work1); double* p_h = mp_work2; const double *p_Ykc_to_use = p_Ykc; std::vector default_Ykc(ne); - if (p_Ykc == NULL) { + if (p_Ykc == nullptr) { for (int i = 0; i < ne; ++i) { default_Ykc[i] = 0.0; } @@ -1119,10 +1119,10 @@ void Thermodynamics::surfaceMassBalance( sum += p_Xw[i]; } - // Use "large" amount of condences phase to simulate infinite surface - double LargeNumber = 100.0; + // Use "large" amount of condensed phase to simulate infinite surface + const double LargeNumber = 100.0; std::vector condensedPhaseElements; - double tol = 1.0e-16; + const double tol = 1.0e-16; for (int i = 0; i < ne; ++i) { p_Xw[i] += LargeNumber*p_Ykc_to_use[i]; sum += LargeNumber*p_Ykc_to_use[i]; diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index 92f7df64..f66f8385 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -886,7 +886,7 @@ class Thermodynamics //: public StateModelUpdateHandler void surfaceMassBalance( const double *const p_Yke, const double *const p_Ykg, const double T, const double P, const double Bg, double &Bc, double &hw, - double *const p_Xs = NULL, const double *const p_Ykc = NULL); + double *const p_Xs = nullptr, const double *const p_Ykc = nullptr); private: