initial work on Radiation1D class

Author: wandadars

include/cantera/oneD/Flow1D.h

@@ -7,11 +7,13 @@
 #define CT_FLOW1D_H
 
 #include "Domain1D.h"
+#include "Radiation1D.h"
 #include "cantera/base/Array.h"
 #include "cantera/base/Solution.h"
 #include "cantera/thermo/ThermoPhase.h"
 #include "cantera/kinetics/Kinetics.h"
 
+
 namespace Cantera
 {
 
@@ -259,12 +261,12 @@ class Flow1D : public Domain1D
 
     //! Return emissivity at left boundary
     double leftEmissivity() const {
-        return m_epsilon_left;
+        return m_radiation->leftEmissivity();
     }
 
     //! Return emissivity at right boundary
     double rightEmissivity() const {
-        return m_epsilon_right;
+        return m_radiation->rightEmissivity();
     }
 
     //! Specify that the the temperature should be held fixed at point `j`.
@@ -461,38 +463,8 @@ class Flow1D : public Domain1D
     //! to be updated are defined.
     virtual void updateProperties(size_t jg, double* x, size_t jmin, size_t jmax);
 
-    /**
-     * Computes the radiative heat loss vector over points jmin to jmax and stores
-     * the data in the qdotRadiation variable.
-     *
-     * The `fit-type` of `polynomial` is uses the model described below.
-     *
-     * The simple radiation model used was established by Liu and Rogg
-     * @cite liu1991. This model considers the radiation of CO2 and H2O.
-     *
-     * This model uses the optically thin limit and the gray-gas approximation to
-     * simply calculate a volume specified heat flux out of the Planck absorption
-     * coefficients, the boundary emissivities and the temperature. Polynomial lines
-     * calculate the species Planck coefficients for H2O and CO2. The data for the
-     * lines are taken from the RADCAL program @cite RADCAL.
-     * The coefficients for the polynomials are taken from
-     * [TNF Workshop](https://tnfworkshop.org/radiation/) material.
-     *
-     *
-     * The `fit-type` of `table` is uses the model described below.
-     *
-     * Spectra for molecules are downloaded with HAPI library from // https://hitran.org/hapi/
-     * [R.V. Kochanov, I.E. Gordon, L.S. Rothman, P. Wcislo, C. Hill, J.S. Wilzewski,
-     * HITRAN Application Programming Interface (HAPI): A comprehensive approach
-     * to working with spectroscopic data, J. Quant. Spectrosc. Radiat. Transfer 177,
-     * 15-30 (2016), https://doi.org/10.1016/j.jqsrt.2016.03.005].
-     *
-     * Planck mean optical path lengths are what are read in from a YAML input file.
-     *
-     *
-     *
-     */
-    void computeRadiation(double* x, size_t jmin, size_t jmax);
+    //! Compute the radiative heat loss at each grid point
+    void computeRadiation(double*, size_t, size_t);
 
     //! @}
 
@@ -918,6 +890,9 @@ class Flow1D : public Domain1D
     //! radiative heat loss.
     double m_epsilon_right = 0.0;
 
+    //! Radiation object used for calculating radiative heat loss
+    std::unique_ptr<Radiation1D> m_radiation;
+
     //! Indices within the ThermoPhase of the radiating species. First index is
     //! for CO2, second is for H2O.
     vector<size_t> m_kRadiating;
@@ -968,16 +943,6 @@ class Flow1D : public Domain1D
     //! radiative heat loss vector
     vector<double> m_qdotRadiation;
 
-    // boundary emissivities for the radiation calculations
-    double m_epsilon_left = 0.0;
-    double m_epsilon_right = 0.0;
-
-    std::map<std::string, int> m_absorptionSpecies; //!< Absorbing species
-    AnyMap m_PMAC; //!< Absorption coefficient data for each species
-
-    // Old radiation variable that can not be deleted for some reason
-    std::vector<size_t> m_kRadiating;
-
     // fixed T and Y values
     //! Fixed values of the temperature at each grid point that are used when solving
     //! with the energy equation disabled.
@@ -1032,5 +997,4 @@ class Flow1D : public Domain1D
 };
 
 }
-
 #endif

include/cantera/oneD/Radiation1D.h

@@ -0,0 +1,105 @@
+//! @file Radiation1D.h
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef RADIATION1D_H
+#define RADIATION1D_H
+
+#include "Domain1D.h"
+#include "cantera/base/Array.h"
+#include "cantera/thermo/ThermoPhase.h"
+
+#include <functional>
+
+namespace Cantera
+{
+
+/**
+ * Computes the radiative heat loss vector over points jmin to jmax and stores
+ * the data in the qdotRadiation variable.
+ *
+ * The `fit-type` of `polynomial` is uses the model described below.
+ *
+ * The simple radiation model used was established by Liu and Rogg
+ * @cite liu1991. This model considers the radiation of CO2 and H2O.
+ *
+ * This model uses the optically thin limit and the gray-gas approximation to
+ * simply calculate a volume specified heat flux out of the Planck absorption
+ * coefficients, the boundary emissivities and the temperature. Polynomial lines
+ * calculate the species Planck coefficients for H2O and CO2. The data for the
+ * lines are taken from the RADCAL program @cite RADCAL.
+ * The coefficients for the polynomials are taken from
+ * [TNF Workshop](https://tnfworkshop.org/radiation/) material.
+ *
+ *
+ * The `fit-type` of `table` is uses the model described below.
+ *
+ * Spectra for molecules are downloaded with HAPI library from // https://hitran.org/hapi/
+ * [R.V. Kochanov, I.E. Gordon, L.S. Rothman, P. Wcislo, C. Hill, J.S. Wilzewski,
+ * HITRAN Application Programming Interface (HAPI): A comprehensive approach
+ *  to working with spectroscopic data, J. Quant. Spectrosc. Radiat. Transfer 177,
+ *   15-30 (2016), https://doi.org/10.1016/j.jqsrt.2016.03.005].
+ *
+ * Planck mean optical path lengths are what are read in from a YAML input file.
+*/
+class Radiation1D {
+public:
+    Radiation1D(ThermoPhase* thermo, double pressure, size_t points,
+                std::function<double(const double*, size_t)> temperatureFunction,
+                std::function<double(const double*, size_t, size_t)> moleFractionFunction);
+
+    // Parse radiation data from YAML input
+    void parseRadiationData();
+
+    // Compute radiative heat loss
+    void computeRadiation(double* x, size_t jmin, size_t jmax, vector<double>& qdotRadiation);
+
+     //! Set the emissivities for the boundary values
+    /*!
+     * Reads the emissivities for the left and right boundary values in the
+     * radiative term and writes them into the variables, which are used for the
+     * calculation.
+     */
+    void setBoundaryEmissivities(double e_left, double e_right);
+
+    //! Return emissivity at left boundary
+    double leftEmissivity() const {
+        return m_epsilon_left;
+    }
+
+    //! Return emissivity at right boundary
+    double rightEmissivity() const {
+        return m_epsilon_right;
+    }
+
+private:
+    ThermoPhase* m_thermo;
+    double m_press;
+    size_t m_points;
+
+    map<string, size_t> m_absorptionSpecies; //!< Absorbing species
+    AnyMap m_PMAC; //!< Absorption coefficient data for each species
+
+    //! Emissivity of the surface to the left of the domain. Used for calculating
+    //! radiative heat loss.
+    double m_epsilon_left = 0.0;
+
+    //! Emissivity of the surface to the right of the domain. Used for calculating
+    //! radiative heat loss.
+    double m_epsilon_right = 0.0;
+
+    // Helper functions
+    double calculatePolynomial(const vector<double>& coefficients, double temperature);
+    double interpolateTable(const vector<double>& temperatures, const vector<double>& data, double temperature);
+
+    //! Lambda function to get temperature at a given point
+    std::function<double(const double*, size_t)> m_T;
+
+    //! Lambda function to get mole fraction at a given point
+    std::function<double(const double*, size_t, size_t)> m_X;
+};
+
+} // namespace Cantera
+
+#endif // RADIATION1D_H

src/oneD/Flow1D.cpp

@@ -20,7 +20,10 @@ namespace Cantera
 
 Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) :
     Domain1D(nsp+c_offset_Y, points),
-    m_nsp(nsp)
+    m_nsp(nsp),
+    m_radiation(make_unique<Radiation1D>(ph, ph->pressure(), points,
+                [this](const double* x, size_t j) {return this->T(x,j);},
+                [this](const double* x, size_t k, size_t j) {return this->X(x,k,j);}))
 {
     m_points = points;
 
@@ -82,70 +85,6 @@ Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) :
         gr.push_back(1.0*ng/m_points);
     }
     setupGrid(m_points, gr.data());
-
-    // Parse radiation data from the YAML file
-    for (auto& name : m_thermo->speciesNames()) {
-        auto& data = m_thermo->species(name)->input;
-        if (data.hasKey("radiation")) {
-            cout << "Radiation data found for species " << name << endl;
-            m_absorptionSpecies.insert({name, m_thermo->speciesIndex(name)});
-            if (data["radiation"].hasKey("fit-type")) {
-                m_PMAC[name]["fit-type"] = data["radiation"]["fit-type"].asString();
-            } else {
-                throw InputFileError("Flow1D::Flow1D", data,
-                "No 'fit-type' entry found for species '{}'", name);
-            }
-
-            // This is the direct tabulation of the optical path length
-            if (data["radiation"]["fit-type"] == "table") {
-                if (data["radiation"].hasKey("temperatures")) {
-                    cout << "Storing temperatures for species " << name << endl;
-                    // Each species may have a specific set of temperatures that are used
-                    m_PMAC[name]["temperatures"] = data["radiation"]["temperatures"].asVector<double>();
-                } else {
-                    throw InputFileError("Flow1D::Flow1D", data,
-                    "No 'temperatures' entry found for species '{}'", name);
-                }
-                if (data["radiation"].hasKey("data")) {
-                    cout << "Storing data for species " << name << endl;
-                    // This data is the Plank mean absorption coefficient
-                    m_PMAC[name]["coefficients"] = data["radiation"]["data"].asVector<double>();
-                } else {
-                    throw InputFileError("Flow1D::Flow1D", data,
-                    "No 'data' entry found for species '{}'", name);
-                }
-            } else if (data["radiation"]["fit-type"] == "polynomial") {
-                cout << "Polynomial fit found for species " << name << endl;
-                if (data["radiation"].hasKey("data")) {
-                    cout << "Storing data for species " << name << endl;
-                    m_PMAC[name]["coefficients"] = data["radiation"]["data"].asVector<double>();
-                } else {
-                    throw InputFileError("Flow1D::Flow1D", data,
-                    "No 'data' entry found for species '{}'", name);
-                }
-            } else {
-                throw InputFileError("Flow1D::Flow1D", data,
-                "Invalid 'fit-type' entry found for species '{}'", name);
-            }
-        }
-    }
-
-    // Polynomial coefficients for CO2 and H2O (backwards compatibility)
-    // Check if "CO2" is already in the map, if not, add the polynomial fit data
-    if (!m_PMAC.hasKey("CO2")) {
-        const std::vector<double> c_CO2 = {18.741, -121.310, 273.500, -194.050, 56.310,
-                                           -5.8169};
-        m_PMAC["CO2"]["fit-type"] = "polynomial";
-        m_PMAC["CO2"]["coefficients"] = c_CO2;
-    }
-
-    // Check if "H2O" is already in the map, if not, add the polynomial fit data
-    if (!m_PMAC.hasKey("H2O")) {
-        const std::vector<double> c_H2O = {-0.23093, -1.12390, 9.41530, -2.99880,
-                                           0.51382, -1.86840e-5};
-        m_PMAC["H2O"]["fit-type"] = "polynomial";
-        m_PMAC["H2O"]["coefficients"] = c_H2O;
-    }
 }
 
 Flow1D::Flow1D(shared_ptr<ThermoPhase> th, size_t nsp, size_t points)
@@ -527,70 +466,7 @@ void Flow1D::updateDiffFluxes(const double* x, size_t j0, size_t j1)
 
 void Flow1D::computeRadiation(double* x, size_t jmin, size_t jmax)
 {
-    // Variable definitions for the Planck absorption coefficient and the
-    // radiation calculation:
-    double k_P_ref = 1.0*OneAtm;
-
-    // Calculation of the two boundary values
-    double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4);
-    double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4);
-
-    double coef = 0.0;
-    for (size_t j = jmin; j < jmax; j++) {
-        // calculation of the mean Planck absorption coefficient
-        double k_P = 0;
-
-        for(const auto& [sp_name, sp_idx] : m_absorptionSpecies) {
-            if (m_PMAC[sp_name]["fit-type"].asString() == "table") {
-                // temperature table interval search
-                int T_index = 0;
-                const int OPL_table_size = m_PMAC[sp_name]["temperatures"].asVector<double>().size();
-                for (int k = 0; k < OPL_table_size; k++) {
-                    if (T(x, j) < m_PMAC[sp_name]["temperatures"].asVector<double>()[k]) {
-                        if (T(x, j) < m_PMAC[sp_name]["temperatures"].asVector<double>()[0]) {
-                            T_index = 0; //lower table limit
-                        }
-                        else {
-                            T_index = k;
-                        }
-                        break;
-                    }
-                    else {
-                        T_index=OPL_table_size-1; //upper table limit
-                    }
-                }
-
-                // absorption coefficient for specie
-                double k_P_specie = 0.0;
-                if ((T_index == 0) || (T_index == OPL_table_size-1)) {
-                    coef=log(1.0/m_PMAC[sp_name]["coefficients"].asVector<double>()[T_index]);
-                }
-                else {
-                    coef=log(1.0/m_PMAC[sp_name]["coefficients"].asVector<double>()[T_index-1])+
-                    (log(1.0/m_PMAC[sp_name]["coefficients"].asVector<double>()[T_index])-log(1.0/m_PMAC[sp_name]["coefficients"].asVector<double>()[T_index-1]))*
-                    (T(x, j)-m_PMAC[sp_name]["temperatures"].asVector<double>()[T_index-1])/(m_PMAC[sp_name]["temperatures"].asVector<double>()[T_index]-m_PMAC[sp_name]["temperatures"].asVector<double>()[T_index-1]);
-                }
-                k_P_specie = exp(coef);
-
-                k_P_specie /= k_P_ref;
-                k_P += m_press * X(x, m_absorptionSpecies[sp_name], j) * k_P_specie;
-            } else if (m_PMAC[sp_name]["fit-type"].asString() == "polynomial") {
-                double k_P_specie = 0.0;
-                for (size_t n = 0; n <= 5; n++) {
-                    k_P_specie += m_PMAC[sp_name]["coefficients"].asVector<double>()[n] * pow(1000 / T(x, j), (double) n);
-                }
-                k_P_specie /= k_P_ref;
-                k_P += m_press * X(x, m_absorptionSpecies[sp_name], j) * k_P_specie;
-            }
-        }
-
-        // Calculation of the radiative heat loss term
-        double radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4)
-                                     - boundary_Rad_left - boundary_Rad_right);
-
-        // set the radiative heat loss vector
-        m_qdotRadiation[j] = radiative_heat_loss;
-    }
+    m_radiation->computeRadiation(x, jmin, jmax, m_qdotRadiation);
 }
 
 void Flow1D::evalContinuity(double* x, double* rsd, int* diag,
@@ -1182,16 +1058,7 @@ bool Flow1D::doElectricField(size_t j) const
 
 void Flow1D::setBoundaryEmissivities(double e_left, double e_right)
 {
-    if (e_left < 0 || e_left > 1) {
-        throw CanteraError("Flow1D::setBoundaryEmissivities",
-            "The left boundary emissivity must be between 0.0 and 1.0!");
-    } else if (e_right < 0 || e_right > 1) {
-        throw CanteraError("Flow1D::setBoundaryEmissivities",
-            "The right boundary emissivity must be between 0.0 and 1.0!");
-    } else {
-        m_epsilon_left = e_left;
-        m_epsilon_right = e_right;
-    }
+    m_radiation->setBoundaryEmissivities(e_left, e_right);
 }
 
 void Flow1D::fixTemperature(size_t j)