Skip to content

Draft: Closes #255. Hold NASA Cp constant when temp is out-of-bounds. #257

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mgoodson-cvd wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Draft: Closes #255. Hold NASA Cp constant when temp is out-of-bounds. #257

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c2ef1d0
Hold NASA Cp constant when temp is out-of-bounds.
mgoodson-cvd Mar 7, 2024
4383174
Update contributors list.
mgoodson-cvd Mar 7, 2024
6869155
Update version to 1.2.0.
mgoodson-cvd Mar 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
)

Expand Down
2 changes: 2 additions & 0 deletions docs/contributors.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@ Pierre Schrooyen <br>
Laurent Soucasse <br>
Alessandro Turchi <br>
[Ruben Di Battista](https://rdb.is) <br>
Matthew Goodson - [Corvid Technologies](https://www.corvidtec.com)
Alireza Mazaheri - NASA LaRC

Have you contributed but don't see your name? Open a PR on Github and add
yourself in this file, we will be happy to add you.
Expand Down
169 changes: 122 additions & 47 deletions src/thermo/NasaDB.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,51 +49,51 @@ class NasaDB : public ThermoDB
NasaDB()
: ThermoDB(298.15, 101325.0), m_ns(0)
{ }

virtual ~NasaDB() {};

virtual bool speciesThermoValidAtT(const size_t i, const double T) const {
assert(i < m_ns);
return (T > m_polynomials[i].minT() && T <= m_polynomials[i].maxT());
}

void cp(
double Th, double Te, double Tr, double Tv, double Tel,
double* const cp, double* const cpt, double* const cpr,
double Th, double Te, double Tr, double Tv, double Tel,
double* const cp, double* const cpt, double* const cpr,
double* const cpv, double* const cpel);

void enthalpy(
double Th, double Te, double Tr, double Tv, double Tel,
double* const h, double* const ht, double* const hr, double* const hv,
double Th, double Te, double Tr, double Tv, double Tel,
double* const h, double* const ht, double* const hr, double* const hv,
double* const hel, double* const hf);

void entropy(
double Th, double Te, double Tr, double Tv, double Tel, double P,
double* const s, double* const st, double* const sr, double* const sv,
double* const s, double* const st, double* const sr, double* const sv,
double* const sel);

void gibbs(
double Th, double Te, double Tr, double Tv, double Tel, double P,
double* const g, double* const gt, double* const gr, double* const gv,
double* const g, double* const gt, double* const gr, double* const gv,
double* const gel);

protected:

void loadAvailableSpecies(std::list<Species>& species_list);

void loadThermodynamicData();

virtual std::string filename() const = 0;

virtual void skipHeader(std::ifstream& is) const = 0;

virtual Species loadSpecies(std::ifstream& is) const = 0;

virtual void loadPolynomials(
std::ifstream& is, std::vector<PolynomialType>& polynomials) = 0;

private:

size_t m_ns;
std::vector<PolynomialType> m_polynomials;
double mp_params[8];
Expand All @@ -102,53 +102,97 @@ class NasaDB : public ThermoDB
template <typename PolynomialType>
void NasaDB<PolynomialType>::cp(
double Th, double Te, double Tr, double Tv, double Tel,
double* const cp, double* const cpt, double* const cpr,
double* const cp, double* const cpt, double* const cpr,
double* const cpv, double* const cpel)
{

if (cp != NULL) {
PolynomialType::computeParams(Th, mp_params, PolynomialType::CP);
for (size_t i = 0; i < m_ns; ++i)
for (size_t i = 0; i < m_ns; ++i) {
// Limit temperature of fit to min and max of polynomial data
double Tfit = std::min(std::max(Th, m_polynomials[i].minT()), m_polynomials[i].maxT());
PolynomialType::computeParams(Tfit, mp_params, PolynomialType::CP);
m_polynomials[i].cp(mp_params, cp[i]);
}
}

if (cpt != NULL) std::fill(cpt, cpt+m_ns, 0.0);
if (cpr != NULL) std::fill(cpr, cpr+m_ns, 0.0);
if (cpv != NULL) std::fill(cpv, cpv+m_ns, 0.0);
if (cpel != NULL) std::fill(cpel, cpel+m_ns, 0.0);
}

template <typename PolynomialType>
void NasaDB<PolynomialType>::enthalpy(
double Th, double Te, double Tr, double Tv, double Tel,
double* const h, double* const ht, double* const hr, double* const hv,
double* const h, double* const ht, double* const hr, double* const hv,
double* const hel, double* const hf)
{
if (h != NULL) {
PolynomialType::computeParams(Th, mp_params, PolynomialType::ENTHALPY);
for (size_t i = 0; i < m_ns; ++i)
for (size_t i = 0; i < m_ns; ++i) {
// Limit temperature of fit to min and max of polynomial data
double const Tmin = m_polynomials[i].minT();
double const Tmax = m_polynomials[i].maxT();
double Tfit = Th;
bool out_of_range = false;
if (Th < Tmin) {
Tfit = Tmin;
out_of_range = true;
} else if (Th > Tmax) {
Tfit = Tmax;
out_of_range = true;
}
PolynomialType::computeParams(Tfit, mp_params, PolynomialType::ENTHALPY);
m_polynomials[i].enthalpy(mp_params, h[i]);
// If outside of range, correct using constant Cp
if (out_of_range) {
h[i] *= Tfit/Th; // Non-dimensionalize by Th
PolynomialType::computeParams(Tfit, mp_params, PolynomialType::CP);
double cp;
m_polynomials[i].cp(mp_params, cp);
h[i] += cp*(Th-Tfit)/Th;
}
}
}

if (ht != NULL) std::fill(ht, ht+m_ns, 0.0);
if (hr != NULL) std::fill(hr, hr+m_ns, 0.0);
if (hv != NULL) std::fill(hv, hv+m_ns, 0.0);
if (hel != NULL) std::fill(hel, hel+m_ns, 0.0);
if (hf != NULL) std::fill(hf, hf+m_ns, 0.0);
if (hf != NULL) std::fill(hf, hf+m_ns, 0.0);
}

template <typename PolynomialType>
void NasaDB<PolynomialType>::entropy(
double Th, double Te, double Tr, double Tv, double Tel, double P,
double* const s, double* const st, double* const sr, double* const sv,
double* const s, double* const st, double* const sr, double* const sv,
double* const sel)
{
if (s != NULL) {
PolynomialType::computeParams(Th, mp_params, PolynomialType::ENTROPY);
for (size_t i = 0; i < m_ns; ++i)
for (size_t i = 0; i < m_ns; ++i) {
// Limit temperature of fit to min and max of polynomial data
double const Tmin = m_polynomials[i].minT();
double const Tmax = m_polynomials[i].maxT();
double Tfit = Th;
bool out_of_range = false;
if (Th < Tmin) {
Tfit = Tmin;
out_of_range = true;
} else if (Th > Tmax) {
Tfit = Tmax;
out_of_range = true;
}
PolynomialType::computeParams(Tfit, mp_params, PolynomialType::ENTROPY);
m_polynomials[i].entropy(mp_params, s[i]);
// If outside of range, correct using constant Cp
if (out_of_range) {
PolynomialType::computeParams(Tfit, mp_params, PolynomialType::CP);
double cp;
m_polynomials[i].cp(mp_params, cp);
s[i] += cp*(log(Th)-log(Tfit));
}
}
}

if (st != NULL) std::fill(st, st+m_ns, 0.0);
if (sr != NULL) std::fill(sr, sr+m_ns, 0.0);
if (sv != NULL) std::fill(sv, sv+m_ns, 0.0);
Expand All @@ -158,15 +202,46 @@ void NasaDB<PolynomialType>::entropy(
template <typename PolynomialType>
void NasaDB<PolynomialType>::gibbs(
double Th, double Te, double Tr, double Tv, double Tel, double P,
double* const g, double* const gt, double* const gr, double* const gv,
double* const g, double* const gt, double* const gr, double* const gv,
double* const gel)
{
if (g != NULL) {
PolynomialType::computeParams(Th, mp_params, PolynomialType::GIBBS);
for (size_t i = 0; i < m_ns; ++i)
m_polynomials[i].gibbs(mp_params, g[i]);
for (size_t i = 0; i < m_ns; ++i) {
// Limit temperature of fit to min and max of polynomial data
double const Tmin = m_polynomials[i].minT();
double const Tmax = m_polynomials[i].maxT();
double Tfit = Th;
bool out_of_range = false;
if (Th < Tmin) {
Tfit = Tmin;
out_of_range = true;
} else if (Th > Tmax) {
Tfit = Tmax;
out_of_range = true;
}
if (out_of_range) {
// Compute G from corrected H and S using constant Cp
PolynomialType::computeParams(Tfit, mp_params, PolynomialType::ENTHALPY);
double h;
m_polynomials[i].enthalpy(mp_params, h);
h *= Tfit/Th; // Non-dimensionalize by Th
double s;
PolynomialType::computeParams(Tfit, mp_params, PolynomialType::ENTROPY);
m_polynomials[i].entropy(mp_params, s);
PolynomialType::computeParams(Tfit, mp_params, PolynomialType::CP);
double cp;
m_polynomials[i].cp(mp_params, cp);
h += cp*(Th-Tfit)/Th;
s += cp*(log(Th)-log(Tfit));
g[i] = h - s;
} else {
// Directly compute G from combined H and S polynomial
PolynomialType::computeParams(Th, mp_params, PolynomialType::GIBBS);
m_polynomials[i].gibbs(mp_params, g[i]);
}
}
}

if (gt != NULL) std::fill(gt, gt+m_ns, 0.0);
if (gr != NULL) std::fill(gr, gr+m_ns, 0.0);
if (gv != NULL) std::fill(gv, gv+m_ns, 0.0);
Expand All @@ -186,16 +261,16 @@ void NasaDB<PolynomialType>::loadAvailableSpecies(
throw FileNotFoundError(db_path)
<< "Could not find thermodynamic database.";
}

// Move to the beginning of the first species
skipHeader(file);

// Load all of the available species until we reach the end of the file
species_list.push_back(loadSpecies(file));
while (species_list.back().name() != "")
species_list.push_back(loadSpecies(file));
species_list.pop_back();

// Close the database file
file.close();
}
Expand All @@ -207,20 +282,20 @@ void NasaDB<PolynomialType>::loadThermodynamicData()
std::string db_path =
Utilities::databaseFileName(filename(), "thermo", ".dat");
std::ifstream file(db_path.c_str());

if (!file.is_open()) {
throw FileNotFoundError(db_path)
<< "Could not find thermodynamic database.";
}

// Move to the beginning of the first species
skipHeader(file);

// Ask the concrete class to fill the polynomials
m_ns = species().size();
m_polynomials.resize(m_ns);
loadPolynomials(file, m_polynomials);

// Close the database file
file.close();
}
Expand All @@ -229,4 +304,4 @@ void NasaDB<PolynomialType>::loadThermodynamicData()
} // namespace Mutation