From d32a72c568da73edbfa6b49d7eba6fdc0bf20693 Mon Sep 17 00:00:00 2001 From: Yue Wang <yue.wang@iws.uni-stuttgart.de> Date: Thu, 23 Mar 2023 13:40:33 +0100 Subject: [PATCH] [bin/util] merge two files in make_co2_table the script generates a directly includable file. --- bin/util/co2table.hh.template | 198 ++++++++++++++++++++++++++++++++ bin/util/co2values.inc.template | 60 ---------- bin/util/make_co2_table.py | 4 +- 3 files changed, 200 insertions(+), 62 deletions(-) create mode 100644 bin/util/co2table.hh.template delete mode 100644 bin/util/co2values.inc.template diff --git a/bin/util/co2table.hh.template b/bin/util/co2table.hh.template new file mode 100644 index 0000000000..5c2c9021a9 --- /dev/null +++ b/bin/util/co2table.hh.template @@ -0,0 +1,198 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program 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 General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup Components + * \brief A generic template for tabulated material laws that depend + * on two parameters. + */ +#ifndef DUMUX_COMPONENTS_CO2_CO2TABLES_HH +#define DUMUX_COMPONENTS_CO2_CO2TABLES_HH + +#include <dune/common/float_cmp.hh> + + /* Tables for CO2 fluid properties calculated according to Span and + * Wagner (1996) and using a web service of the National Institute + * of Standards and Techlology of the U.S. Department of Commerce: + * https://webbook.nist.gov/ + * + * THIS AN AUTO-GENERATED FILE! DO NOT EDIT IT! + * + ******************************************************************** + + In case you are using this the data generated with this script + please cite the following publications: + + P.J. Linstrom and W.G. Mallard, Eds., + NIST Chemistry WebBook, NIST Standard Reference Database Number 69, + National Institute of Standards and Technology, Gaithersburg MD, 20899, + https://doi.org/10.18434/T4D303, (retrieved $DATE). + + Span, Roland, and Wolfgang Wagner. + "A new equation of state for carbon dioxide covering + the fluid region from the tripleâ€point temperature + to 1100 K at pressures up to 800 MPa." + Journal of physical and chemical reference data 25.6 (1996): 1509-1596. + https://doi.org/10.1063/1.555991 + + ******************************************************************** + * + * Generated using: + * + * ./make_co2_table.py -t1 $MIN_TEMP -t2 $MAX_TEMP -nt $NUM_TEMP_SAMPLES -p1 $MIN_PRESS -p2 $MAX_PRESS -np $NUM_PRESS_SAMPLES + */ + + +namespace Dumux::GeneratedCO2Tables { + +/*! + * \ingroup Components + * \brief A generic template for tabulated material laws that depend + * on two parameters. + */ +template <class Traits> +class TabulatedProperties +{ + using Scalar = typename Traits::Scalar; + + static constexpr auto numTempSteps = Traits::numTempSteps; + static constexpr auto numPressSteps = Traits::numPressSteps; + +public: + TabulatedProperties() = default; + + constexpr Scalar minTemp() const { return Traits::minTemp; } + constexpr Scalar maxTemp() const { return Traits::maxTemp; } + constexpr Scalar minPress() const { return Traits::minPress; } + constexpr Scalar maxPress() const { return Traits::maxPress; } + + constexpr bool applies(Scalar temperature, Scalar pressure) const + { + return minTemp() <= temperature && temperature <= maxTemp() && + minPress() <= pressure && pressure <= maxPress(); + } + + constexpr Scalar at(Scalar temperature, Scalar pressure) const + { + if (!applies(temperature, pressure)) + { + if (temperature<minTemp()) temperature = minTemp(); + else if (temperature>maxTemp()) temperature = maxTemp(); + + if (pressure<minPress()) pressure = minPress(); + else if (pressure>maxPress()) pressure = maxPress(); + } + + const int i = findTempIdx_(temperature); + const int j = findPressIdx_(pressure); + + const Scalar tempAtI = temperatureAt_(i); + const Scalar tempAtI1 = temperatureAt_(i + 1); + const Scalar pressAtI = pressureAt_(j); + const Scalar pressAtI1 = pressureAt_(j + 1); + + const Scalar alpha = (temperature - tempAtI)/(tempAtI1 - tempAtI); + const Scalar beta = (pressure - pressAtI)/(pressAtI1 - pressAtI); + + // bi-linear interpolation + const Scalar lowresValue = + (1-alpha)*(1-beta)*val(i, j) + + (1-alpha)*( beta)*val(i, j + 1) + + ( alpha)*(1-beta)*val(i + 1, j) + + ( alpha)*( beta)*val(i + 1, j + 1); + + // return the weighted sum of the low- and high-resolution values + return lowresValue; + } + + constexpr Scalar val(int i, int j) const + { return Traits::vals[i][j]; } + +private: + constexpr int findTempIdx_(Scalar temperature) const + { + if (Dune::FloatCmp::eq<Scalar>(temperature, maxTemp())) + return numTempSteps - 2; + + const int result = static_cast<int>((temperature - minTemp())/(maxTemp() - minTemp())*(numTempSteps - 1)); + + using std::clamp; + return clamp(result, 0, numTempSteps - 2); + } + + constexpr int findPressIdx_(Scalar pressure) const + { + if (Dune::FloatCmp::eq<Scalar>(pressure, maxPress())) + return numPressSteps - 2; + + const int result = static_cast<int>((pressure - minPress())/(maxPress() - minPress())*(numPressSteps - 1)); + + using std::clamp; + return clamp(result, 0, numPressSteps - 2); + } + + constexpr Scalar temperatureAt_(int i) const + { return i*(maxTemp() - minTemp())/(numTempSteps - 1) + minTemp(); } + + constexpr Scalar pressureAt_(int j) const + { return j*(maxPress() - minPress())/(numPressSteps - 1) + minPress(); } +}; + +struct TabulatedDensityTraits +{ + using Scalar = double; + static constexpr std::string_view name = "density"; + static constexpr int numTempSteps = $NUM_TEMP_SAMPLES; + static constexpr Scalar minTemp = $MIN_TEMP; + static constexpr Scalar maxTemp = $MAX_TEMP; + static constexpr int numPressSteps = $NUM_PRESS_SAMPLES; + static constexpr Scalar minPress = $MIN_PRESS; + static constexpr Scalar maxPress = $MAX_PRESS; + static constexpr Scalar vals[numTempSteps][numPressSteps] = { + $DENSITY_VALS + }; +}; + +struct TabulatedEnthalpyTraits +{ + using Scalar = double; + static constexpr std::string_view name = "enthalpy"; + static constexpr int numTempSteps = $NUM_TEMP_SAMPLES; + static constexpr Scalar minTemp = $MIN_TEMP; + static constexpr Scalar maxTemp = $MAX_TEMP; + static constexpr int numPressSteps = $NUM_PRESS_SAMPLES; + static constexpr Scalar minPress = $MIN_PRESS; + static constexpr Scalar maxPress = $MAX_PRESS; + static constexpr Scalar vals[numTempSteps][numPressSteps] = { + $ENTHALPY_VALS + }; +}; + +using TabulatedDensity = TabulatedProperties<TabulatedDensityTraits>; +using TabulatedEnthalpy = TabulatedProperties<TabulatedEnthalpyTraits>; + +// this class collects all the tabulated quantities in one convenient place +struct CO2Tables +{ + static constexpr inline TabulatedEnthalpy tabulatedEnthalpy = {}; + static constexpr inline TabulatedDensity tabulatedDensity = {}; +}; +} // end namespace Dumux::GeneratedCO2Tables + +#endif diff --git a/bin/util/co2values.inc.template b/bin/util/co2values.inc.template deleted file mode 100644 index 9f16b1ab2e..0000000000 --- a/bin/util/co2values.inc.template +++ /dev/null @@ -1,60 +0,0 @@ -/* Tables for CO2 fluid properties calculated according to Span and - * Wagner (1996) and using a web service of the National Institute - * of Standards and Techlology of the U.S. Department of Commerce: - * https://webbook.nist.gov/ - * - * THIS AN AUTO-GENERATED FILE! DO NOT EDIT IT! - * - ******************************************************************** - - In case you are using this the data generated with this script - please cite the following publications: - - P.J. Linstrom and W.G. Mallard, Eds., - NIST Chemistry WebBook, NIST Standard Reference Database Number 69, - National Institute of Standards and Technology, Gaithersburg MD, 20899, - https://doi.org/10.18434/T4D303, (retrieved $DATE). - - Span, Roland, and Wolfgang Wagner. - "A new equation of state for carbon dioxide covering - the fluid region from the tripleâ€point temperature - to 1100 K at pressures up to 800 MPa." - Journal of physical and chemical reference data 25.6 (1996): 1509-1596. - https://doi.org/10.1063/1.555991 - - ******************************************************************** - * - * Generated using: - * - * ./make_co2_table.py -t1 $MIN_TEMP -t2 $MAX_TEMP -nt $NUM_TEMP_SAMPLES -p1 $MIN_PRESS -p2 $MAX_PRESS -np $NUM_PRESS_SAMPLES - */ - -struct TabulatedDensityTraits -{ - using Scalar = double; - static constexpr std::string_view name = "density"; - static constexpr int numTempSteps = $NUM_TEMP_SAMPLES; - static constexpr Scalar minTemp = $MIN_TEMP; - static constexpr Scalar maxTemp = $MAX_TEMP; - static constexpr int numPressSteps = $NUM_PRESS_SAMPLES; - static constexpr Scalar minPress = $MIN_PRESS; - static constexpr Scalar maxPress = $MAX_PRESS; - static constexpr Scalar vals[numTempSteps][numPressSteps] = { - $DENSITY_VALS - }; -}; - -struct TabulatedEnthalpyTraits -{ - using Scalar = double; - static constexpr std::string_view name = "enthalpy"; - static constexpr int numTempSteps = $NUM_TEMP_SAMPLES; - static constexpr Scalar minTemp = $MIN_TEMP; - static constexpr Scalar maxTemp = $MAX_TEMP; - static constexpr int numPressSteps = $NUM_PRESS_SAMPLES; - static constexpr Scalar minPress = $MIN_PRESS; - static constexpr Scalar maxPress = $MAX_PRESS; - static constexpr Scalar vals[numTempSteps][numPressSteps] = { - $ENTHALPY_VALS - }; -}; diff --git a/bin/util/make_co2_table.py b/bin/util/make_co2_table.py index e00ecb2a4b..d436b39ae3 100755 --- a/bin/util/make_co2_table.py +++ b/bin/util/make_co2_table.py @@ -152,7 +152,7 @@ DENSITY_DATA = ",\n".join(DENSITY_DATA) ENTHALPY_DATA = ",\n".join(ENTHALPY_DATA) # write the table by filling the gaps in the template -with open("co2values.inc.template", "r") as templateFile: +with open("co2table.hh.template", "r") as templateFile: template = Template(templateFile.read()) replacements = { @@ -167,5 +167,5 @@ replacements = { "DATE": date.today().strftime("%B %d, %Y"), } -with open("co2values.inc", "w") as tables: +with open("co2table.hh", "w") as tables: tables.write(template.substitute(replacements)) -- GitLab