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