8#include <compile_time_options.h>
12#include <deal.II/base/array_view.h>
13#include <deal.II/base/config.h>
14#include <deal.II/base/exceptions.h>
15#include <deal.II/base/parameter_acceptor.h>
18#include "eos_Interface.h"
38 enum class TableType : EOS_INTEGER {
60 Interface(
const std::vector<std::tuple<EOS_INTEGER, TableType>> &tables)
62 n_tables_ = tables.size();
64 std::transform(std::begin(tables),
66 std::back_inserter(material_ids_),
67 [](
const auto &it) {
return std::get<0>(it); });
69 std::transform(std::begin(tables),
71 std::back_inserter(table_types_),
73 return static_cast<EOS_INTEGER
>(std::get<1>(it));
76 table_handles_.resize(n_tables_);
80 EOS_INTEGER error_code;
81 eos_CreateTables(&n_tables_,
84 table_handles_.data(),
86 check_tables(
"eos_CreateTables");
90 for (EOS_INTEGER i = 0; i < n_tables_; i++) {
93 &table_handles_[i], &EOS_SMOOTH, EOS_NullPtr, &error_code);
94 check_error_code(error_code,
"eos_SetOption", i);
99 eos_LoadTables(&n_tables_, table_handles_.data(), &error_code);
100 check_tables(
"eos_LoadTables");
106 ~Interface() noexcept
108 EOS_INTEGER error_code;
109 eos_DestroyTables(&n_tables_, table_handles_.data(), &error_code);
116 inline DEAL_II_ALWAYS_INLINE
void
117 interpolate_values(
const EOS_INTEGER &index,
118 const dealii::ArrayView<EOS_REAL> &F,
119 const dealii::ArrayView<EOS_REAL> &dFx,
120 const dealii::ArrayView<EOS_REAL> &dFy,
121 const dealii::ArrayView<const EOS_REAL> &X,
122 const dealii::ArrayView<const EOS_REAL> &Y)
124 Assert(index >= 0 && index < n_tables_,
125 dealii::ExcMessage(
"Table index out of range"));
127 EOS_INTEGER n_queries = F.size();
129 Assert(dFx.size() == n_queries && dFy.size() == n_queries &&
130 X.size() == n_queries && Y.size() == n_queries,
131 dealii::ExcMessage(
"vector sizes do not match"));
133 EOS_INTEGER error_code;
134 eos_Interpolate(&table_handles_[index],
136 const_cast<EOS_REAL *
>(X.data()),
137 const_cast<EOS_REAL *
>(Y.data()),
148 std::vector<EOS_INTEGER> material_ids_;
149 std::vector<EOS_INTEGER> table_types_;
150 EOS_INTEGER n_tables_;
152 mutable std::vector<EOS_INTEGER> table_handles_;
158 void check_error_code(
159 EOS_INTEGER error_code,
160 const std::string &routine,
161 EOS_INTEGER i = std::numeric_limits<EOS_INTEGER>::max())
const
163 if (error_code != EOS_OK) {
164 std::array<EOS_CHAR, EOS_MaxErrMsgLen> error_message;
165 eos_GetErrorMessage(&error_code, error_message.data());
167 std::stringstream exception_body;
168 exception_body <<
"Error: " << routine;
169 if (i != std::numeric_limits<EOS_INTEGER>::max())
170 exception_body <<
" (table " << i <<
")";
171 exception_body <<
": " << error_code <<
" - " << error_message.data()
173 AssertThrow(
false, dealii::ExcMessage(exception_body.str()));
177 void check_tables(
const std::string &routine)
const
179 for (EOS_INTEGER i = 0; i < n_tables_; i++) {
180 EOS_INTEGER table_error_code = EOS_OK;
181 eos_GetErrorCode(&table_handles_[i], &table_error_code);
182 if (table_error_code != EOS_OK) {
183 std::array<EOS_CHAR, EOS_MaxErrMsgLen> error_message;
184 eos_GetErrorMessage(&table_error_code, error_message.data());
186 std::stringstream exception_body;
187 exception_body <<
"Error: " << routine <<
" (table " << i
188 <<
"): " << table_error_code <<
" - "
189 << error_message.data() << std::flush;
191 AssertThrow(
false, dealii::ExcMessage(exception_body.str()));
199 namespace EquationOfStateLibrary
216 Sesame(
const std::string &subsection)
221 "material id", material_id_,
"The Sesame Material ID");
225 const auto set_up_database = [&]() {
226 const std::vector<std::tuple<EOS_INTEGER, eospac::TableType>> tables{
227 {material_id_, eospac::TableType::p_rho_e},
228 {material_id_, eospac::TableType::e_rho_p},
231 eospac_interface_ = std::make_unique<eospac::Interface>(tables);
234 this->parse_parameters_call_back.connect(set_up_database);
238 double pressure(
double rho,
double e)
const final
240 EOS_INTEGER index = 0;
242 double p, p_drho, p_de;
243 const double rho_scaled = rho / 1.0e3;
244 const double e_scaled = e / 1.0e6;
246 eospac_interface_->interpolate_values(
248 dealii::ArrayView<double>(&p, 1),
249 dealii::ArrayView<double>(&p_drho, 1),
250 dealii::ArrayView<double>(&p_de, 1),
251 dealii::ArrayView<const double>(&rho_scaled, 1),
252 dealii::ArrayView<const double>(&e_scaled, 1));
258 void pressure(
const dealii::ArrayView<double> &p,
259 const dealii::ArrayView<double> &rho,
260 const dealii::ArrayView<double> &e)
const final
262 Assert(p.size() == rho.size() && rho.size() == e.size(),
263 dealii::ExcMessage(
"vectors have different size"));
265 EOS_INTEGER index = 0;
268 thread_local static std::vector<double> p_drho;
269 thread_local static std::vector<double> p_de;
270 p_drho.resize(p.size());
271 p_de.resize(p.size());
274 std::transform(std::begin(rho),
277 [](
double rho) {
return rho / 1.0e3; });
280 std::transform(std::begin(e),
283 [](
auto e) {
return e / 1.0e6; });
285 eospac_interface_->interpolate_values(index,
287 dealii::ArrayView<double>(p_drho),
288 dealii::ArrayView<double>(p_de),
293 std::transform(std::begin(p),
296 [](
auto it) {
return it * 1.0e9; });
302 EOS_INTEGER index = 1;
304 double e, e_drho, e_dp;
305 const double rho_scaled = rho / 1.0e3;
306 const double p_scaled = p / 1.0e9;
308 eospac_interface_->interpolate_values(
310 dealii::ArrayView<double>(&e, 1),
311 dealii::ArrayView<double>(&e_drho, 1),
312 dealii::ArrayView<double>(&e_dp, 1),
313 dealii::ArrayView<const double>(&rho_scaled, 1),
314 dealii::ArrayView<const double>(&p_scaled, 1));
322 const dealii::ArrayView<double> &rho,
323 const dealii::ArrayView<double> &p)
const final
325 Assert(e.size() == rho.size() && rho.size() == p.size(),
326 dealii::ExcMessage(
"vectors have different size"));
328 EOS_INTEGER index = 1;
331 thread_local static std::vector<double> e_drho;
332 thread_local static std::vector<double> e_dp;
333 e_drho.resize(e.size());
334 e_dp.resize(e.size());
337 std::transform(std::begin(rho),
340 [](
double rho) {
return rho / 1.0e3; });
343 std::transform(std::begin(p),
346 [](
auto it) {
return it / 1.0e9; });
348 eospac_interface_->interpolate_values(index,
350 dealii::ArrayView<double>(e_drho),
351 dealii::ArrayView<double>(e_dp),
356 std::transform(std::begin(e),
359 [](
auto e) {
return e * 1.0e6; });
369 void temperature(
const dealii::ArrayView<double> & ,
370 const dealii::ArrayView<double> & ,
371 const dealii::ArrayView<double> & )
const final
383 const dealii::ArrayView<double> & ,
384 const dealii::ArrayView<double> & )
const final
390 EOS_INTEGER material_id_;
391 std::unique_ptr<eospac::Interface> eospac_interface_;
402 "ryujin has to be configured with eospac support in order to use "
403 "the Sesame EOS database";
407 AssertThrow(
false, dealii::ExcMessage(
message));
413 AssertThrow(
false, dealii::ExcMessage(
message));
419 AssertThrow(
false, dealii::ExcMessage(
message));
425 AssertThrow(
false, dealii::ExcMessage(
message));
EquationOfState(const std::string &name, const std::string &subsection)
bool prefer_vector_interface_
double pressure(double, double) const final
double speed_of_sound(double, double) const final
double temperature(double, double) const final
Sesame(const std::string &subsection)
double specific_internal_energy(double, double) const final
static constexpr auto message