Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/EoS/ThermallyPerfectGas.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Thermally perfect gas equation of state
9 : : \details This file declares functions for the thermally perfect gas equation
10 : : of state for the compressible flow equations.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef ThermallyPerfectGas_h
14 : : #define ThermallyPerfectGas_h
15 : :
16 : : #include "Data.hpp"
17 : :
18 : : namespace inciter {
19 : :
20 [ + - ]: 48 : class ThermallyPerfectGas {
21 : :
22 : : private:
23 : : tk::real m_R;
24 : : std::vector< std::vector< tk::real > > m_cp_coeff{3, std::vector< tk::real >(8)};
25 : : std::vector< tk::real > m_t_range{std::vector< tk::real >(4)};
26 : : tk::real m_dH_ref;
27 : :
28 : 35373742 : void get_t_range( tk::real &temp_poly,
29 : : std::size_t &t_rng_idx ) const
30 : : // *************************************************************************
31 : : //! \brief Check what temperature range the given temperature is in. If it
32 : : //! exceeds the bounds, reset the temp to the bounds.
33 : : //! \param[in] temp Given temperature to be checked for range
34 : : //! \return Index of temperature range the given temperature is in
35 : : // *************************************************************************
36 : : {
37 : : // First, bounds check
38 [ - + ]: 35373742 : if (temp_poly < m_t_range[0]) {
39 : 0 : t_rng_idx = 0;
40 : 0 : temp_poly = m_t_range[0];
41 [ + - ]: 35373742 : } else if (temp_poly > m_t_range.back()) {
42 : 0 : t_rng_idx = m_t_range.size() - 1;
43 : 0 : temp_poly = m_t_range.back();
44 : : } else {
45 : : // Valid bounds
46 [ + - ]: 52114523 : for (std::size_t k = 0; k < m_t_range.size() - 1; k++) {
47 [ + - ][ + + ]: 52114523 : if (temp_poly >= m_t_range[k] && temp_poly <= m_t_range[k+1]) {
48 : 35373742 : t_rng_idx = k;
49 : 35373742 : break;
50 : : }
51 : : }
52 : : }
53 : 35373742 : }
54 : :
55 : 15240329 : tk::real calc_h(tk::real temp) const
56 : : // *************************************************************************
57 : : //! Calculate dimensionless enthalpy according to the NASA-9 polynomial
58 : : //! \param[in] temp temperature at which to calculate enthalpy
59 : : //! \return dimensionless enthalpy, h / (R * T)
60 : : // *************************************************************************
61 : : {
62 : : // Identify what temperature range this falls in. If it falls outside the
63 : : // bounds, some corrections must be applied.
64 : 15240329 : std::size_t t_rng_idx(0); // Reference the correct polynomial.
65 : 15240329 : tk::real temp_poly = temp; // temp to use in polynomial expression
66 : 15240329 : get_t_range(temp_poly, t_rng_idx); // Bounds check performed inside
67 : :
68 : : // h = h_poly(T) + h_ref
69 [ - + ]: 15240329 : tk::real h = -m_cp_coeff[t_rng_idx][0] * std::pow(temp_poly, -2) +
70 : 15240329 : m_cp_coeff[t_rng_idx][1] * std::log(temp_poly) / temp_poly +
71 : 15240329 : m_cp_coeff[t_rng_idx][2] +
72 : 15240329 : m_cp_coeff[t_rng_idx][3] * temp_poly / 2. +
73 : 15240329 : m_cp_coeff[t_rng_idx][4] * std::pow(temp_poly, 2) / 3. +
74 : 15240329 : m_cp_coeff[t_rng_idx][5] * std::pow(temp_poly, 3) / 4. +
75 : 15240329 : m_cp_coeff[t_rng_idx][6] * std::pow(temp_poly, 4) / 5. +
76 : 15240329 : m_cp_coeff[t_rng_idx][7] / temp_poly;
77 : :
78 : : // If bounds exceeded, temp_poly will be different than temp. Apply correction.
79 [ - + ]: 15240329 : if (std::abs(temp_poly - temp) > std::numeric_limits< tk::real >::epsilon()) {
80 : 0 : tk::real cp_star = calc_cp(temp_poly);
81 : 0 : h = h * temp_poly / temp + (temp - temp_poly) / temp * cp_star;
82 : : }
83 : :
84 : 15240329 : return h;
85 : : }
86 : :
87 : 20133413 : tk::real calc_cp(tk::real temp) const
88 : : // *************************************************************************
89 : : //! Calculate dimensionless specific heat according to the NASA-9 polynomial
90 : : //! \param[in] temp temperature at which to calculate specific heat
91 : : //! \return dimensionless enthalpy, c_p / R
92 : : // *************************************************************************
93 : : {
94 : : // Identify what temperature range this falls in. If it falls outside the
95 : : // bounds, some corrections must be applied.
96 : 20133413 : std::size_t t_rng_idx(0); // Reference the correct polynomial.
97 : 20133413 : tk::real temp_poly = temp; // temp to use in polynomial expression
98 : 20133413 : get_t_range(temp_poly, t_rng_idx); // Bounds check performed inside
99 : :
100 : 20133413 : tk::real cp = m_cp_coeff[t_rng_idx][0] * std::pow(temp_poly, -2) +
101 : 20133413 : m_cp_coeff[t_rng_idx][1] / temp_poly +
102 : 20133413 : m_cp_coeff[t_rng_idx][2] +
103 : 20133413 : m_cp_coeff[t_rng_idx][3] * temp_poly +
104 : 20133413 : m_cp_coeff[t_rng_idx][4] * std::pow(temp_poly, 2) +
105 : 20133413 : m_cp_coeff[t_rng_idx][5] * std::pow(temp_poly, 3) +
106 : 20133413 : m_cp_coeff[t_rng_idx][6] * std::pow(temp_poly, 4);
107 : :
108 : 20133413 : return cp;
109 : : }
110 : :
111 : : public:
112 : : //! Default constructor
113 : : ThermallyPerfectGas() = default;
114 : :
115 : : //! Constructor
116 : : ThermallyPerfectGas(
117 : : tk::real R,
118 : : std::vector< std::vector< tk::real > > cp_coeff,
119 : : std::vector< tk::real > t_range,
120 : : tk::real dH_ref);
121 : :
122 : : //! Set rho0 EOS parameter. No-op.
123 : : void setRho0(tk::real) {}
124 : :
125 : : //! Calculate density from the material pressure and temperature
126 : : [[noreturn]] tk::real density( tk::real pr,
127 : : tk::real temp ) const;
128 : :
129 : : //! Calculate pressure from the material density, momentum and total energy
130 : : [[noreturn]] tk::real pressure( tk::real rho,
131 : : tk::real u,
132 : : tk::real v,
133 : : tk::real w,
134 : : tk::real rhoE,
135 : : tk::real alpha=1.0,
136 : : std::size_t imat=0,
137 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad={{}}) const;
138 : :
139 : : //! Calculate cold-compression component of pressure (no-op)
140 : : tk::real pressure_coldcompr(
141 : : tk::real,
142 : : tk::real ) const
143 : : { return 0.0; }
144 : :
145 : : //! \brief Calculate the Cauchy stress tensor from the material density,
146 : : //! momentum, and total energy
147 : : std::array< std::array< tk::real, 3 >, 3 >
148 : : CauchyStress(
149 : : tk::real,
150 : : tk::real,
151 : : tk::real,
152 : : tk::real,
153 : : tk::real,
154 : : tk::real,
155 : : std::size_t,
156 : : const std::array< std::array< tk::real, 3 >, 3 >& adefgrad={{}} ) const;
157 : :
158 : : //! Calculate speed of sound from the material density and material pressure
159 : : [[noreturn]] tk::real soundspeed( tk::real rho,
160 : : tk::real pr,
161 : : tk::real alpha=1.0,
162 : : std::size_t imat=0,
163 : : const std::array< std::array< tk::real, 3 >, 3 >& adefgrad={{}},
164 : : const std::array< tk::real, 3 >& asigman={{}} ) const;
165 : :
166 : : //! Calculate speed of shear waves
167 : : tk::real shearspeed(
168 : : tk::real,
169 : : tk::real,
170 : : std::size_t ) const { return 0.0; }
171 : :
172 : : //! \brief Calculate material specific total energy from the material
173 : : //! density, momentum and material pressure
174 : : [[noreturn]] tk::real totalenergy( tk::real rho,
175 : : tk::real u,
176 : : tk::real v,
177 : : tk::real w,
178 : : tk::real pr,
179 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad={{}} ) const;
180 : :
181 : : //! \brief Calculate material temperature from the material density, and
182 : : //! material specific total energy
183 : : [[noreturn]] tk::real temperature( tk::real rho,
184 : : tk::real u,
185 : : tk::real v,
186 : : tk::real w,
187 : : tk::real rhoE,
188 : : tk::real alpha=1.0,
189 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad={{}} ) const;
190 : :
191 : : //! Compute the minimum allowed pressure
192 : : tk::real min_eff_pressure(
193 : : tk::real min,
194 : : tk::real,
195 : : tk::real ) const
196 : : { return min; }
197 : :
198 : : //! Compute the reference density
199 : : tk::real refDensity() const { return density(refPressure(), 300.0); }
200 : :
201 : : //! Compute the reference pressure
202 : : tk::real refPressure() const { return 1.0e5; }
203 : :
204 : : //! Return initial density
205 : 0 : tk::real rho0() const { return density(1.0e5, 300.0); }
206 : :
207 : : //! Return gas constant (species specific)
208 : : tk::real gas_constant() const { return m_R; }
209 : :
210 : : //! Return species internal energy
211 : : tk::real internalenergy(tk::real temp) const;
212 : :
213 : : //! Return species specific heat (constant volume)
214 : : tk::real cv(tk::real temp) const;
215 : :
216 : : /** @name Charm++ pack/unpack serializer member functions */
217 : : ///@{
218 : : //! \brief Pack/Unpack serialize member function
219 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
220 : : void pup( PUP::er &p ) /*override*/ {
221 : : p | m_R;
222 : : p | m_cp_coeff;
223 : : p | m_t_range;
224 : : p | m_dH_ref;
225 : : }
226 : : //! \brief Pack/Unpack serialize operator|
227 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
228 : : //! \param[in,out] i ThermallyPerfectGas object reference
229 : : friend void operator|( PUP::er& p, ThermallyPerfectGas& i ) { i.pup(p); }
230 : : //@}
231 : : };
232 : :
233 : : } //inciter::
234 : :
235 : : #endif // ThermallyPerfectGas_h
|
