diff --git a/.gitignore b/.gitignore index c7921345c33..e9df3e32dc6 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ __pycache__ *.exe.manifest build/ test/work/ +test/data/cross include/cantera/base/config.h include/cantera/base/config.h.build include/cantera/base/system.h.gch diff --git a/include/cantera/kinetics/Photolysis.h b/include/cantera/kinetics/Photolysis.h index 2ce58e6feda..3d70edd4f31 100644 --- a/include/cantera/kinetics/Photolysis.h +++ b/include/cantera/kinetics/Photolysis.h @@ -64,11 +64,10 @@ class PhotolysisBase : public ReactionRate { void setParameters(AnyMap const& node, UnitStack const& rate_units) override; - void setRateParameters(const AnyValue& rate, vector const& branches); - - void loadCrossSectionVulcan(vector files, string const& branch); - - void loadCrossSectionKinetics7(vector files, string const& branch); + //! Set the rate parameters for each branch + //! @param rate Rate coefficient data + //! @param branch_map Map of branch names to branch indices + void setRateParameters(const AnyValue& rate, map const& branch_map); void getParameters(AnyMap& node) const override; @@ -126,23 +125,14 @@ class PhotolysisRate : public PhotolysisBase { } double evalFromStruct(PhotolysisData const& data) { - if (m_crossSection.empty()) { - return 0.; - } - double wmin = m_temp_wave_grid[m_ntemp]; double wmax = m_temp_wave_grid.back(); - if (wmin > data.wavelength.front()) { - throw CanteraError("PhotolysisRate::evalFromStruct", - "Wavelength out of range: {} nm < {} nm", - wmin, data.wavelength.front()); - } - - if (wmax < data.wavelength.back()) { - throw CanteraError("PhotolysisRate::evalFromStruct", - "Wavelength out of range: {} nm > {} nm", - wmax, data.wavelength.back()); + if (m_crossSection.empty() || + wmin > data.wavelength.back() || + wmax < data.wavelength.front()) + { + return 0.; } int iwmin = locate(data.wavelength.data(), wmin, data.wavelength.size()); @@ -158,8 +148,9 @@ class PhotolysisRate : public PhotolysisBase { len, 2, m_branch.size()); double total_rate = 0.0; - for (auto& [name, stoich] : m_net_products) - stoich = 0.0; + for (auto const& branch : m_branch) + for (auto const& [name, stoich] : branch) + m_net_products[name] = 0.; for (int i = iwmin; i < iwmax; i++) { coord[1] = data.wavelength[i+1]; @@ -169,8 +160,9 @@ class PhotolysisRate : public PhotolysisBase { for (size_t n = 0; n < m_branch.size(); n++) { double rate = 0.5 * (data.wavelength[i+1] - data.wavelength[i]) * (cross1[n] * data.actinicFlux[i] + cross2[n] * data.actinicFlux[i+1]); - for (auto const& [name, stoich] : m_branch[n]) - m_net_products[name] += rate * stoich; + for (auto const& [name, stoich] : m_branch[n]) { + m_net_products.at(name) += rate * stoich; + } total_rate += rate; cross1[n] = cross2[n]; } @@ -190,6 +182,23 @@ class PhotolysisRate : public PhotolysisBase { Composition m_net_products; }; +/** + * Read the cross-section data from VULCAN format files + * + * @param files Vector of filenames. + * There are two files for each photolysis reaction. The first one is for + * cross-section data and the second one for the branch ratios. + */ +pair, vector> +load_xsection_vulcan(vector const& files, vector const& branches); + +/** + * Read the cross-section data from KINETICS7 format files + */ + +pair, vector> +load_xsection_kinetics7(vector const& files, vector const& branches); + } #endif // CT_PHOTOLYSIS_H diff --git a/src/kinetics/Photolysis.cpp b/src/kinetics/Photolysis.cpp index 75cfe0febc6..2d6be0827b3 100644 --- a/src/kinetics/Photolysis.cpp +++ b/src/kinetics/Photolysis.cpp @@ -76,7 +76,8 @@ bool PhotolysisData::check() const } PhotolysisBase::PhotolysisBase( - vector const& temp, vector const& wavelength, + vector const& temp, + vector const& wavelength, vector const& branches, vector const& xsection): m_crossSection(xsection) @@ -111,7 +112,8 @@ PhotolysisBase::PhotolysisBase(AnyMap const& node, UnitStack const& rate_units) setParameters(node, rate_units); } -void PhotolysisBase::setRateParameters(const AnyValue& rate, vector const& branches) +void PhotolysisBase::setRateParameters(AnyValue const& rate, + map const& branch_map) { if (rate.hasKey("resolution")) { double resolution = rate["resolution"].asDouble(); @@ -129,13 +131,13 @@ void PhotolysisBase::setRateParameters(const AnyValue& rate, vector cons } } else { for (auto const& [branch, scale] : rate["scale"].as()) { - auto b = std::find(branches.begin(), branches.end(), branch); - if (b == branches.end()) { + auto it = branch_map.find(branch); + if (it == branch_map.end()) { throw CanteraError("PhotolysisBase::setRateParameters", "Branch '{}' not found", branch); } - scales[b - branches.begin()] = scale.asDouble(); + scales[it->second] = scale.asDouble(); } } } @@ -143,20 +145,18 @@ void PhotolysisBase::setRateParameters(const AnyValue& rate, vector cons void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_units) { - vector branches; + map branch_map; + pair, vector> result; vector temperature; - vector wavelength; - vector xsection; ReactionRate::setParameters(node, rate_units); if (node.hasKey("branches")) { for (auto const& branch : node["branches"].asVector()) { - branches.push_back(branch["name"].asString()); + branch_map[branch["name"].asString()] = m_branch.size(); m_branch.push_back(parseCompString(branch["product"].asString())); } } else { - branches.push_back("all"); vector tokens; tokenizeString(node["equation"].asString(), tokens); m_branch.push_back(parseCompString(tokens[0] + ":1")); @@ -165,7 +165,6 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni if (node.hasKey("cross-section")) { for (auto const& data: node["cross-section"].asVector()) { auto format = data["format"].asString(); - auto branch = data.hasKey("branch") ? data["branch"].asString() : "all"; auto temp = data["temperature-range"].asVector(2, 2); if (temp[0] >= temp[1]) { throw CanteraError("PhotolysisBase::setParameters", @@ -186,15 +185,15 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni if (format == "YAML") { for (auto const& entry: data["data"].asVector>()) { - wavelength.push_back(entry[0]); - xsection.push_back(entry[1]); + result.first.push_back(entry[0]); + result.second.push_back(entry[1]); } } else if (format == "VULCAN") { auto files = data["filenames"].asVector(); - loadCrossSectionVulcan(files, branch); + result = load_xsection_vulcan(files, m_branch); } else if (format == "KINETICS7") { auto files = data["filenames"].asVector(); - loadCrossSectionKinetics7(files, branch); + result = load_xsection_kinetics7(files, m_branch); } else { throw CanteraError("PhotolysisBase::setParameters", "Invalid cross-section format '{}'.", format); @@ -202,6 +201,9 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni } } + vector &wavelength(result.first); + vector &xsection(result.second); + m_ntemp = temperature.size(); m_nwave = wavelength.size(); m_temp_wave_grid.resize(m_ntemp + m_nwave); @@ -219,7 +221,7 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni m_crossSection.insert(m_crossSection.end(), xsection.begin(), xsection.end()); if (node.hasKey("rate-constant")) { - setRateParameters(node["rate-constant"], branches); + setRateParameters(node["rate-constant"], branch_map); } m_valid = true; @@ -262,12 +264,4 @@ void PhotolysisBase::validate(string const& equation, Kinetics const& kin) } } -void __attribute__((weak)) PhotolysisBase::loadCrossSectionVulcan(vector files, - string const& branch) -{} - -void __attribute__((weak)) PhotolysisBase::loadCrossSectionKinetics7(vector files, - string const& branch) -{} - } diff --git a/src/kinetics/load_xsection_kinetics7.cpp b/src/kinetics/load_xsection_kinetics7.cpp new file mode 100644 index 00000000000..62821c4759b --- /dev/null +++ b/src/kinetics/load_xsection_kinetics7.cpp @@ -0,0 +1,120 @@ +#include +#include + +#include "cantera/base/stringUtils.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/Photolysis.h" + +namespace Cantera +{ + +pair, vector> +load_xsection_kinetics7(vector const& files, vector const& branches) +{ + if (files.size() != 1) { + throw CanteraError("load_xsection_kinetics7", + "Only one file can be loaded for Kinetics7 format."); + } + + auto const& filename = files[0]; + + FILE* file = fopen(filename.c_str(), "r"); + + if (!file) { + throw CanteraError("load_xsection_kinetics7", + "Could not open file '{}'", filename); + } + + vector wavelength; + vector xsection; + + int nbranch = branches.size(); + int min_is = 9999, max_ie = 0; + + char *line = NULL; + size_t len = 0; + ssize_t read; + + // Read each line from the file + while ((read = getline(&line, &len, file)) != -1) { + // Skip empty lines or lines containing only whitespace + if (line[0] == '\n') continue; + + char equation[61]; + int is, ie, nwave; + float temp; + + // read header + int num = sscanf(line, "%60c%4d%4d%4d%6f", equation, &is, &ie, &nwave, &temp); + min_is = std::min(min_is, is); + max_ie = std::max(max_ie, ie); + + if (num != 5) { + throw CanteraError("PhotolysisBase::loadCrossSectionKinetics7", + "Header format from file '{}' is wrong.", filename); + } + + // initialize wavelength and xsection for the first time + if (wavelength.size() == 0) { + wavelength.resize(nwave); + xsection.resize(nwave * nbranch); + } + + // read content + int ncols = 7; + int nrows = ceil(1. * nwave / ncols); + + auto product = parseCompString(equation); + + auto it = std::find(branches.begin(), branches.end(), product); + + if (it == branches.end()) { + // skip this section + for (int i = 0; i < nrows; i++) + getline(&line, &len, file); + } else { + for (int i = 0; i < nrows; i++) { + getline(&line, &len, file); + + for (int j = 0; j < ncols; j++) { + float wave, cross; + int num = sscanf(line + 17*j, "%7f%10f", &wave, &cross); + if (num != 2) { + throw CanteraError("PhotolysisBase::loadCrossSectionKinetics7", + "Cross-section format from file '{}' is wrong.", filename); + } + int b = it - branches.begin(); + int k = i * ncols + j; + + if (k >= nwave) break; + wavelength[k] = wave / 10.; // Angstrom to nm + xsection[k * nbranch + b] = cross * 10.; // cm^2 / Angstrom to cm^2 / nm + } + } + } + } + + // remove unused wavelength and xsection + wavelength = vector(wavelength.begin() + min_is - 1, + wavelength.begin() + max_ie); + + xsection = vector(xsection.begin() + (min_is - 1) * nbranch, + xsection.begin() + max_ie * nbranch); + + // A -> A is the total cross section in kinetics7 format + // need to subtract the other branches + + for (size_t i = 0; i < wavelength.size(); i++) { + for (int j = 1; j < nbranch; j++) { + xsection[i * nbranch] -= xsection[i * nbranch + j]; + } + xsection[i * nbranch] = std::max(xsection[i * nbranch], 0.); + } + + free(line); + fclose(file); + + return {std::move(wavelength), std::move(xsection)}; +} + +} // namespace Cantera diff --git a/src/kinetics/load_xsection_vulcan.cpp b/src/kinetics/load_xsection_vulcan.cpp new file mode 100644 index 00000000000..40e2e4b9869 --- /dev/null +++ b/src/kinetics/load_xsection_vulcan.cpp @@ -0,0 +1,24 @@ +#include +#include + +#include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/Photolysis.h" + +namespace Cantera +{ + +pair, vector> +load_xsection_vulcan(vector const& files, vector const& branches) +{ + if (files.size() != 2) { + throw CanteraError("load_xsection_Vulcan", + "Only two files can be loaded for Vulcan format."); + } + + vector wavelength; + vector xsection; + + return {std::move(wavelength), std::move(xsection)}; +} + +} // namespace Cantera diff --git a/test/data/ch4_photolysis.yaml b/test/data/ch4_photolysis.yaml index 148fc858b97..47198cbd28c 100644 --- a/test/data/ch4_photolysis.yaml +++ b/test/data/ch4_photolysis.yaml @@ -69,9 +69,9 @@ reactions: - name: b5 # 0.2 product: "CH:1 H2:1 H:1" cross-section: - - format: VULCAN + - format: KINETICS7 temperature-range: [0., 300.] - filenames: [CH4_cross.csv, CH4_branch.csv] + filenames: [../data/cross/CH4.dat2] rate-constants: resolution: 0.1 scale: @@ -83,6 +83,6 @@ reactions: - format: YAML temperature-range: [0., 300.] data: - - [200., 1.e-4] - - [250., 2.e-4] - - [300., 3.e-4] + - [20., 1.e-18] + - [100., 2.e-18] + - [180., 3.e-18] diff --git a/test/photochem/PhotochemTitan.cpp b/test/photochem/PhotochemTitan.cpp index 335435b2fbe..0a708581043 100644 --- a/test/photochem/PhotochemTitan.cpp +++ b/test/photochem/PhotochemTitan.cpp @@ -25,7 +25,7 @@ class PhotochemTitan: public testing::Test { vector actinic_flux(10); for (int i = 0; i < 10; i++) { - wavelength[i] = 200.0 + i * 10.0; + wavelength[i] = 20.0 + i * 20.0; actinic_flux[i] = 1.0; } @@ -51,8 +51,8 @@ TEST_F(PhotochemTitan, check_fwd_rate_constants) { kin->getFwdRateConstants(kfwd.data()); - ASSERT_NEAR(kfwd[0], 0., 1.0e-8); - ASSERT_NEAR(kfwd[1], 0.0171, 1.0e-8); + ASSERT_NEAR(kfwd[0], 3.06820e-14, 1.0e-18); + ASSERT_NEAR(kfwd[1], 3.2e-16, 1.0e-18); } } // namespace Cantera