Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/EoS/EoS.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 Equation of state class
9 : : \details This file defines functions for equations of state for the
10 : : compressible flow equations.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef EoS_h
14 : : #define EoS_h
15 : :
16 : : #include <cmath>
17 : : #include "Data.hpp"
18 : : #include "Inciter/InputDeck/InputDeck.hpp"
19 : :
20 : : namespace inciter {
21 : :
22 : : extern ctr::InputDeck g_inputdeck;
23 : :
24 : : using ncomp_t = kw::ncomp::info::expect::type;
25 : :
26 : : //! Get a property for a material
27 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
28 : : //! \tparam Prop Tag of property required
29 : : //! \param[in] system Equation system index
30 : : //! \param[in] imat Material-id who's property is required. Default is 0, so
31 : : //! that for the single-material system, this argument can be left unspecified
32 : : //! by the calling code
33 : : //! \return Material ratio of specific heats (gamma)
34 : : template< class Eq, class Prop >
35 : : tk::real
36 : : getmatprop( ncomp_t system, std::size_t imat=0 ) {
37 [ - - ][ + + ]: 207972281 : const auto& matprop = g_inputdeck.get< tag::param, Eq, tag::material >()[ system ];
[ + + ][ + + ]
[ - - ][ - - ]
38 : : const auto& map = g_inputdeck.get< tag::param, Eq, tag::matidxmap >();
39 [ + + ][ - + ]: 445527745 : auto meos = map.template get< tag::eosidx >()[ imat ];
[ + + ]
40 : 472590368 : auto midx = map.template get< tag::matidx >()[ imat ];
41 [ - - ][ + + ]: 223619831 : return matprop[ meos ].template get< Prop >()[ midx ];
[ + + ][ + + ]
[ - - ][ - - ]
42 : : }
43 : :
44 : : //! Get the ratio of specific heats (gamma) for a material
45 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
46 : : //! \param[in] system Equation system index
47 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
48 : : //! for the single-material system, this argument can be left unspecified by
49 : : //! the calling code
50 : : //! \return Material ratio of specific heats (gamma)
51 : : template< class Eq >
52 : : tk::real gamma( ncomp_t system, std::size_t imat=0 )
53 : : {
54 : : return getmatprop< Eq, tag::gamma >(system, imat);
55 : : }
56 : :
57 : : //! Get the specific heat at constant volume (cv) for a material
58 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
59 : : //! \param[in] system Equation system index
60 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
61 : : //! for the single-material system, this argument can be left unspecified by
62 : : //! the calling code
63 : : //! \return Material specific heat at constant volume (cv)
64 : : template< class Eq >
65 : : tk::real cv( ncomp_t system, std::size_t imat=0 )
66 : : {
67 : : return getmatprop< Eq, tag::cv >(system, imat);
68 : : }
69 : :
70 : : //! Get the stiffness parameter (pstiff) for a material
71 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
72 : : //! \param[in] system Equation system index
73 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
74 : : //! for the single-material system, this argument can be left unspecified by
75 : : //! the calling code
76 : : //! \return Material stiffness parameter (pstiff)
77 : : template< class Eq >
78 : : tk::real pstiff( ncomp_t system, std::size_t imat=0 )
79 : : {
80 : : return getmatprop< Eq, tag::pstiff >(system, imat);
81 : : }
82 : :
83 : : //! Get the thermal conductivity (k) for a material
84 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
85 : : //! \param[in] system Equation system index
86 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
87 : : //! for the single-material system, this argument can be left unspecified by
88 : : //! the calling code
89 : : //! \return Material thermal conductivity (k)
90 : : template< class Eq >
91 : : tk::real k( ncomp_t system, std::size_t imat=0 )
92 : : {
93 : : return getmatprop< Eq, tag::k >(system, imat);
94 : : }
95 : :
96 : : //! Get the dynamic viscosity (mu) for a material
97 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
98 : : //! \param[in] system Equation system index
99 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
100 : : //! for the single-material system, this argument can be left unspecified by
101 : : //! the calling code
102 : : //! \return Material dynamic viscosity (mu)
103 : : template< class Eq >
104 : : tk::real mu( ncomp_t system, std::size_t imat=0 )
105 : : {
106 : : return getmatprop< Eq, tag::mu >(system, imat);
107 : : }
108 : :
109 : : //! \brief Calculate density from the material pressure and temperature using
110 : : //! the stiffened-gas equation of state
111 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
112 : : //! \param[in] system Equation system index
113 : : //! \param[in] pr Material pressure
114 : : //! \param[in] temp Material temperature
115 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
116 : : //! for the single-material system, this argument can be left unspecified by
117 : : //! the calling code
118 : : //! \return Material density calculated using the stiffened-gas EoS
119 : : template< class Eq >
120 : 2028900 : tk::real eos_density( ncomp_t system,
121 : : tk::real pr,
122 : : tk::real temp,
123 : : std::size_t imat=0 )
124 : : {
125 : : // query input deck to get gamma, p_c, cv
126 : : auto g = gamma< Eq >(system, imat);
127 : : auto p_c = pstiff< Eq >(system, imat);
128 : : auto c_v = cv< Eq >(system, imat);
129 : :
130 : 2028900 : tk::real rho = (pr + p_c) / ((g-1.0) * c_v * temp);
131 : 2028900 : return rho;
132 : : }
133 : :
134 : : //! \brief Calculate pressure from the material density, momentum and total
135 : : //! energy using the stiffened-gas equation of state
136 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
137 : : //! \param[in] system Equation system index
138 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
139 : : //! \param[in] u X-velocity
140 : : //! \param[in] v Y-velocity
141 : : //! \param[in] w Z-velocity
142 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
143 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for the
144 : : //! single-material system, this argument can be left unspecified by the
145 : : //! calling code
146 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
147 : : //! for the single-material system, this argument can be left unspecified by
148 : : //! the calling code
149 : : //! \return Material partial pressure (alpha_k * p_k) calculated using the
150 : : //! stiffened-gas EoS
151 : : #pragma omp declare simd
152 : : template< class Eq >
153 [ - + ]: 185896409 : tk::real eos_pressure( ncomp_t system,
154 : : tk::real arho,
155 : : tk::real u,
156 : : tk::real v,
157 : : tk::real w,
158 : : tk::real arhoE,
159 : : tk::real alpha=1.0,
160 : : std::size_t imat=0 )
161 : : {
162 : : // query input deck to get gamma, p_c
163 : : auto g = gamma< Eq >(system, imat);
164 : : auto p_c = pstiff< Eq >(system, imat);
165 : :
166 : 185896409 : tk::real partpressure = (arhoE - 0.5 * arho * (u*u + v*v + w*w) - alpha*p_c)
167 [ - + ]: 185896409 : * (g-1.0) - alpha*p_c;
168 : :
169 : : // check partial pressure divergence
170 [ - + ]: 185896409 : if (!std::isfinite(partpressure)) {
171 : : std::cout << "Material-id: " << imat << std::endl;
172 : : std::cout << "Volume-fraction: " << alpha << std::endl;
173 : : std::cout << "Partial density: " << arho << std::endl;
174 : : std::cout << "Total energy: " << arhoE << std::endl;
175 : : std::cout << "Velocity: " << u << ", " << v << ", " << w
176 : : << std::endl;
177 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf partial pressure: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
178 : : + std::to_string(partpressure) + ", material volume fraction: " +
179 : : std::to_string(alpha));
180 : : }
181 : :
182 : 185896409 : return partpressure;
183 : : }
184 : :
185 : : //! Calculate speed of sound from the material density and material pressure
186 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
187 : : //! \param[in] system Equation system index
188 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
189 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
190 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for the
191 : : //! single-material system, this argument can be left unspecified by the
192 : : //! calling code
193 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
194 : : //! for the single-material system, this argument can be left unspecified by
195 : : //! the calling code
196 : : //! \return Material speed of sound using the stiffened-gas EoS
197 : : template< class Eq >
198 [ + + ]: 175343137 : tk::real eos_soundspeed( ncomp_t system,
199 : : tk::real arho, tk::real apr,
200 : : tk::real alpha=1.0, std::size_t imat=0 )
201 : : {
202 : : // query input deck to get gamma, p_c
203 : : auto g = gamma< Eq >(system, imat);
204 : : auto p_c = pstiff< Eq >(system, imat);
205 : :
206 [ + + ]: 175343137 : auto p_eff = std::max( 1.0e-15, apr+(alpha*p_c) );
207 : :
208 [ - + ]: 175343137 : tk::real a = std::sqrt( g * p_eff / arho );
209 : :
210 : : // check sound speed divergence
211 [ - + ]: 175343137 : if (!std::isfinite(a)) {
212 : : std::cout << "Material-id: " << imat << std::endl;
213 : : std::cout << "Volume-fraction: " << alpha << std::endl;
214 : : std::cout << "Partial density: " << arho << std::endl;
215 : : std::cout << "Partial pressure: " << apr << std::endl;
216 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
217 : : + std::to_string(a) + ", material volume fraction: " +
218 : : std::to_string(alpha));
219 : : }
220 : :
221 : 175343137 : return a;
222 : : }
223 : :
224 : : //! \brief Calculate material specific total energy from the material density,
225 : : //! momentum and material pressure
226 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
227 : : //! \param[in] system Equation system index
228 : : //! \param[in] rho Material density
229 : : //! \param[in] u X-velocity
230 : : //! \param[in] v Y-velocity
231 : : //! \param[in] w Z-velocity
232 : : //! \param[in] pr Material pressure
233 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
234 : : //! for the single-material system, this argument can be left unspecified by
235 : : //! the calling code
236 : : //! \return Material specific total energy using the stiffened-gas EoS
237 : : template< class Eq >
238 : 12345150 : tk::real eos_totalenergy( ncomp_t system,
239 : : tk::real rho,
240 : : tk::real u,
241 : : tk::real v,
242 : : tk::real w,
243 : : tk::real pr,
244 : : std::size_t imat=0 )
245 : : {
246 : : // query input deck to get gamma, p_c
247 : : auto g = gamma< Eq >(system, imat);
248 : : auto p_c = pstiff< Eq >(system, imat);
249 : :
250 : 12345150 : tk::real rhoE = (pr + p_c) / (g-1.0) + 0.5 * rho * (u*u + v*v + w*w) + p_c;
251 : 12345150 : return rhoE;
252 : : }
253 : :
254 : : //! \brief Calculate material temperature from the material density, and
255 : : //! material specific total energy
256 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
257 : : //! \param[in] system Equation system index
258 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
259 : : //! \param[in] u X-velocity
260 : : //! \param[in] v Y-velocity
261 : : //! \param[in] w Z-velocity
262 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
263 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for the
264 : : //! single-material system, this argument can be left unspecified by the
265 : : //! calling code
266 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
267 : : //! for the single-material system, this argument can be left unspecified by
268 : : //! the calling code
269 : : //! \return Material temperature using the stiffened-gas EoS
270 : : template< class Eq >
271 : 2229900 : tk::real eos_temperature( ncomp_t system,
272 : : tk::real arho,
273 : : tk::real u,
274 : : tk::real v,
275 : : tk::real w,
276 : : tk::real arhoE,
277 : : tk::real alpha=1.0,
278 : : std::size_t imat=0 )
279 : : {
280 : : // query input deck to get p_c, cv
281 : : auto c_v = cv< Eq >(system, imat);
282 : : auto p_c = pstiff< Eq >(system, imat);
283 : :
284 : 2229900 : tk::real t = (arhoE - 0.5 * arho * (u*u + v*v + w*w) - alpha*p_c) / (arho*c_v);
285 : 2229900 : return t;
286 : : }
287 : :
288 : : //! Constrain material partial pressure (alpha_k * p_k)
289 : : //! \tparam Eq Equation type to operate on, e.g., tag::compflow, tag::multimat
290 : : //! \param[in] system Equation system index
291 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
292 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for the
293 : : //! single-material system, this argument can be left unspecified by the
294 : : //! calling code
295 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
296 : : //! for the single-material system, this argument can be left unspecified by
297 : : //! the calling code
298 : : //! \return Constrained material partial pressure (alpha_k * p_k)
299 : : template< class Eq >
300 [ + + ]: 55814890 : tk::real constrain_pressure( ncomp_t system,
301 : : tk::real apr,
302 : : tk::real alpha=1.0,
303 : : std::size_t imat=0 )
304 : : {
305 : : // query input deck to get p_c
306 : : auto p_c = pstiff< Eq >(system, imat);
307 : :
308 [ + + ]: 55814890 : return std::max(apr, alpha*(-p_c+1e-12));
309 : : }
310 : :
311 : : } //inciter::
312 : :
313 : : #endif // EoS_h
|
