Skip to content

Commit

Permalink
Add KINETICS7 crossdisk (#4)
Browse files Browse the repository at this point in the history
- Add kinetics7 crossdisk
- Check net branch ratio correct

TODO:
- [ ] Connect net branch ratio to Kinetics
  • Loading branch information
chengcli authored Jan 4, 2024
1 parent f18c09a commit 5c42e0f
Show file tree
Hide file tree
Showing 7 changed files with 203 additions and 55 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
55 changes: 32 additions & 23 deletions include/cantera/kinetics/Photolysis.h
Original file line number Diff line number Diff line change
Expand Up @@ -64,11 +64,10 @@ class PhotolysisBase : public ReactionRate {

void setParameters(AnyMap const& node, UnitStack const& rate_units) override;

void setRateParameters(const AnyValue& rate, vector<string> const& branches);

void loadCrossSectionVulcan(vector<string> files, string const& branch);

void loadCrossSectionKinetics7(vector<string> 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<string, int> const& branch_map);

void getParameters(AnyMap& node) const override;

Expand Down Expand Up @@ -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());
Expand All @@ -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];
Expand All @@ -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];
}
Expand All @@ -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<double>, vector<double>>
load_xsection_vulcan(vector<string> const& files, vector<Composition> const& branches);

/**
* Read the cross-section data from KINETICS7 format files
*/

pair<vector<double>, vector<double>>
load_xsection_kinetics7(vector<string> const& files, vector<Composition> const& branches);

}

#endif // CT_PHOTOLYSIS_H
42 changes: 18 additions & 24 deletions src/kinetics/Photolysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,8 @@ bool PhotolysisData::check() const
}

PhotolysisBase::PhotolysisBase(
vector<double> const& temp, vector<double> const& wavelength,
vector<double> const& temp,
vector<double> const& wavelength,
vector<std::string> const& branches,
vector<double> const& xsection):
m_crossSection(xsection)
Expand Down Expand Up @@ -111,7 +112,8 @@ PhotolysisBase::PhotolysisBase(AnyMap const& node, UnitStack const& rate_units)
setParameters(node, rate_units);
}

void PhotolysisBase::setRateParameters(const AnyValue& rate, vector<string> const& branches)
void PhotolysisBase::setRateParameters(AnyValue const& rate,
map<string, int> const& branch_map)
{
if (rate.hasKey("resolution")) {
double resolution = rate["resolution"].asDouble();
Expand All @@ -129,34 +131,32 @@ void PhotolysisBase::setRateParameters(const AnyValue& rate, vector<string> cons
}
} else {
for (auto const& [branch, scale] : rate["scale"].as<AnyMap>()) {
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();
}
}
}
}

void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_units)
{
vector<string> branches;
map<string, int> branch_map;
pair<vector<double>, vector<double>> result;
vector<double> temperature;
vector<double> wavelength;
vector<double> xsection;

ReactionRate::setParameters(node, rate_units);

if (node.hasKey("branches")) {
for (auto const& branch : node["branches"].asVector<AnyMap>()) {
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<string> tokens;
tokenizeString(node["equation"].asString(), tokens);
m_branch.push_back(parseCompString(tokens[0] + ":1"));
Expand All @@ -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<AnyMap>()) {
auto format = data["format"].asString();
auto branch = data.hasKey("branch") ? data["branch"].asString() : "all";
auto temp = data["temperature-range"].asVector<double>(2, 2);
if (temp[0] >= temp[1]) {
throw CanteraError("PhotolysisBase::setParameters",
Expand All @@ -186,22 +185,25 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni

if (format == "YAML") {
for (auto const& entry: data["data"].asVector<vector<double>>()) {
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<string>();
loadCrossSectionVulcan(files, branch);
result = load_xsection_vulcan(files, m_branch);
} else if (format == "KINETICS7") {
auto files = data["filenames"].asVector<string>();
loadCrossSectionKinetics7(files, branch);
result = load_xsection_kinetics7(files, m_branch);
} else {
throw CanteraError("PhotolysisBase::setParameters",
"Invalid cross-section format '{}'.", format);
}
}
}

vector<double> &wavelength(result.first);
vector<double> &xsection(result.second);

m_ntemp = temperature.size();
m_nwave = wavelength.size();
m_temp_wave_grid.resize(m_ntemp + m_nwave);
Expand All @@ -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;
Expand Down Expand Up @@ -262,12 +264,4 @@ void PhotolysisBase::validate(string const& equation, Kinetics const& kin)
}
}

void __attribute__((weak)) PhotolysisBase::loadCrossSectionVulcan(vector<string> files,
string const& branch)
{}

void __attribute__((weak)) PhotolysisBase::loadCrossSectionKinetics7(vector<string> files,
string const& branch)
{}

}
120 changes: 120 additions & 0 deletions src/kinetics/load_xsection_kinetics7.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
#include <cmath>
#include <stdio.h>

#include "cantera/base/stringUtils.h"
#include "cantera/kinetics/Reaction.h"
#include "cantera/kinetics/Photolysis.h"

namespace Cantera
{

pair<vector<double>, vector<double>>
load_xsection_kinetics7(vector<string> const& files, vector<Composition> 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<double> wavelength;
vector<double> 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<double>(wavelength.begin() + min_is - 1,
wavelength.begin() + max_ie);

xsection = vector<double>(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
24 changes: 24 additions & 0 deletions src/kinetics/load_xsection_vulcan.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
#include <cmath>
#include <stdio.h>

#include "cantera/kinetics/Reaction.h"
#include "cantera/kinetics/Photolysis.h"

namespace Cantera
{

pair<vector<double>, vector<double>>
load_xsection_vulcan(vector<string> const& files, vector<Composition> const& branches)
{
if (files.size() != 2) {
throw CanteraError("load_xsection_Vulcan",
"Only two files can be loaded for Vulcan format.");
}

vector<double> wavelength;
vector<double> xsection;

return {std::move(wavelength), std::move(xsection)};
}

} // namespace Cantera
10 changes: 5 additions & 5 deletions test/data/ch4_photolysis.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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]
Loading

0 comments on commit 5c42e0f

Please sign in to comment.