Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "TabulatedFluidProperties.h"
11 : #include "MooseUtils.h"
12 : #include "Conversion.h"
13 : #include "KDTree.h"
14 : #include "BidimensionalInterpolation.h"
15 : #include "NewtonInversion.h"
16 :
17 : // C++ includes
18 : #include <fstream>
19 : #include <ctime>
20 : #include <cmath>
21 : #include <regex>
22 :
23 : InputParameters
24 597 : TabulatedFluidProperties::validParams()
25 : {
26 597 : InputParameters params = SinglePhaseFluidProperties::validParams();
27 597 : params.addClassDescription(
28 : "Single phase fluid properties computed using bi-dimensional interpolation of tabulated "
29 : "values.");
30 :
31 : // Which interpolations to create
32 1194 : params.addParam<bool>("create_pT_interpolations",
33 1194 : true,
34 : "Whether to load (from file) or create (from a fluid property object) "
35 : "properties interpolations from pressure and temperature");
36 1194 : params.addParam<bool>("create_ve_interpolations",
37 1194 : false,
38 : "Whether to load (from file) or create (from a fluid property object) "
39 : "properties interpolations from pressure and temperature");
40 :
41 : // Input / output
42 1194 : params.addParam<UserObjectName>("fp", "The name of the FluidProperties UserObject");
43 1194 : params.addParam<FileName>("fluid_property_file",
44 : "Name of the csv file containing the tabulated fluid property data.");
45 1194 : params.addParam<FileName>(
46 : "fluid_property_ve_file",
47 : "Name of the csv file containing the tabulated (v,e) fluid property data.");
48 1194 : params.addParam<FileName>("fluid_property_output_file",
49 : "Name of the CSV file which can be output with the tabulation. This "
50 : "file can then be read as a 'fluid_property_file'");
51 1194 : params.addParam<FileName>(
52 : "fluid_property_ve_output_file",
53 : "Name of the CSV file which can be output with the (v,e) tabulation. This "
54 : "file can then be read as a 'fluid_property_ve_file'");
55 1194 : params.addDeprecatedParam<bool>(
56 : "save_file",
57 : "Whether to save the csv fluid properties file",
58 : "This parameter is no longer required. Whether to save a CSV tabulation file is controlled "
59 : "by specifying the 'fluid_property_output_file' parameter");
60 :
61 : // Data source on a per-property basis
62 : MultiMooseEnum properties(
63 : "density enthalpy internal_energy viscosity k c cv cp entropy pressure temperature",
64 597 : "density enthalpy internal_energy viscosity");
65 1194 : params.addParam<MultiMooseEnum>("interpolated_properties",
66 : properties,
67 : "Properties to interpolate if no data file is provided");
68 :
69 : // (p,T) grid parameters
70 1194 : params.addRangeCheckedParam<Real>(
71 : "temperature_min", 300, "temperature_min > 0", "Minimum temperature for tabulated data.");
72 1194 : params.addParam<Real>("temperature_max", 500, "Maximum temperature for tabulated data.");
73 1791 : params.addRangeCheckedParam<Real>(
74 1194 : "pressure_min", 1e5, "pressure_min > 0", "Minimum pressure for tabulated data.");
75 1194 : params.addParam<Real>("pressure_max", 50.0e6, "Maximum pressure for tabulated data.");
76 1791 : params.addRangeCheckedParam<unsigned int>(
77 1194 : "num_T", 100, "num_T > 0", "Number of points to divide temperature range.");
78 1791 : params.addRangeCheckedParam<unsigned int>(
79 1194 : "num_p", 100, "num_p > 0", "Number of points to divide pressure range.");
80 :
81 : // (v,e) grid parameters
82 1194 : params.addParam<Real>("e_min", "Minimum specific internal energy for tabulated data.");
83 1194 : params.addParam<Real>("e_max", "Maximum specific internal energy for tabulated data.");
84 1194 : params.addRangeCheckedParam<Real>(
85 : "v_min", "v_min > 0", "Minimum specific volume for tabulated data.");
86 1194 : params.addRangeCheckedParam<Real>(
87 : "v_max", "v_max > 0", "Maximum specific volume for tabulated data.");
88 1194 : params.addParam<bool>("construct_pT_from_ve",
89 1194 : false,
90 : "If the lookup table (p, T) as functions of (v, e) should be constructed.");
91 1194 : params.addParam<bool>("construct_pT_from_vh",
92 1194 : false,
93 : "If the lookup table (p, T) as functions of (v, h) should be constructed.");
94 1791 : params.addRangeCheckedParam<unsigned int>(
95 : "num_v",
96 1194 : 100,
97 : "num_v > 0",
98 : "Number of points to divide specific volume range for (v,e) lookups.");
99 1791 : params.addRangeCheckedParam<unsigned int>("num_e",
100 1194 : 100,
101 : "num_e > 0",
102 : "Number of points to divide specific internal energy "
103 : "range for (v,e) lookups.");
104 1194 : params.addParam<bool>(
105 : "use_log_grid_v",
106 1194 : false,
107 : "Option to use a base-10 logarithmically-spaced grid for specific volume instead of a "
108 : "linearly-spaced grid.");
109 1194 : params.addParam<bool>(
110 : "use_log_grid_e",
111 1194 : false,
112 : "Option to use a base-10 logarithmically-spaced grid for specific internal energy instead "
113 : "of a linearly-spaced grid.");
114 1194 : params.addParam<bool>(
115 : "use_log_grid_h",
116 1194 : false,
117 : "Option to use a base-10 logarithmically-spaced grid for specific enthalpy instead "
118 : "of a linearly-spaced grid.");
119 :
120 : // Out of bounds behavior
121 1194 : params.addDeprecatedParam<bool>(
122 : "error_on_out_of_bounds",
123 : "Whether pressure or temperature from tabulation exceeding user-specified bounds leads to "
124 : "an error.",
125 : "This parameter has been replaced by the 'out_of_bounds_behavior' parameter which offers "
126 : "more flexibility. The option to error is called 'throw' in that parameter.");
127 : // NOTE: this enum must remain the same as OOBBehavior in the header
128 1194 : MooseEnum OOBBehavior("ignore throw declare_invalid warn_invalid set_to_closest_bound", "throw");
129 1194 : params.addParam<MooseEnum>("out_of_bounds_behavior",
130 : OOBBehavior,
131 : "Property evaluation behavior when evaluated outside the "
132 : "user-specified or tabulation-specified bounds");
133 :
134 : // This is generally a bad idea. However, several properties have not been tabulated so several
135 : // tests are relying on the original fp object to provide the value (for example for the
136 : // vaporPressure())
137 1194 : params.addParam<bool>(
138 1194 : "allow_fp_and_tabulation", false, "Whether to allow the two sources of data concurrently");
139 :
140 1194 : params.addParamNamesToGroup("fluid_property_file fluid_property_ve_file "
141 : "fluid_property_output_file fluid_property_ve_output_file",
142 : "Tabulation file read/write");
143 1194 : params.addParamNamesToGroup("construct_pT_from_ve construct_pT_from_vh",
144 : "Variable set conversion");
145 1194 : params.addParamNamesToGroup("temperature_min temperature_max pressure_min pressure_max e_min "
146 : "e_max v_min v_max error_on_out_of_bounds out_of_bounds_behavior",
147 : "Tabulation and interpolation bounds");
148 1194 : params.addParamNamesToGroup(
149 : "num_T num_p num_v num_e use_log_grid_v use_log_grid_e use_log_grid_h",
150 : "Tabulation and interpolation discretization");
151 :
152 597 : return params;
153 597 : }
154 :
155 326 : TabulatedFluidProperties::TabulatedFluidProperties(const InputParameters & parameters)
156 : : SinglePhaseFluidProperties(parameters),
157 652 : _file_name_in(isParamValid("fluid_property_file") ? getParam<FileName>("fluid_property_file")
158 : : ""),
159 354 : _file_name_ve_in(
160 652 : isParamValid("fluid_property_ve_file") ? getParam<FileName>("fluid_property_ve_file") : ""),
161 676 : _file_name_out(isParamValid("fluid_property_output_file")
162 326 : ? getParam<FileName>("fluid_property_output_file")
163 : : ""),
164 652 : _file_name_ve_out(isParamValid("fluid_property_ve_output_file")
165 326 : ? getParam<FileName>("fluid_property_ve_output_file")
166 : : ""),
167 1004 : _save_file(isParamValid("save_file") ? getParam<bool>("save_file")
168 926 : : isParamValid("fluid_property_output_file")),
169 652 : _create_direct_pT_interpolations(getParam<bool>("create_pT_interpolations")),
170 652 : _create_direct_ve_interpolations(getParam<bool>("create_ve_interpolations")),
171 652 : _temperature_min(getParam<Real>("temperature_min")),
172 652 : _temperature_max(getParam<Real>("temperature_max")),
173 652 : _pressure_min(getParam<Real>("pressure_min")),
174 652 : _pressure_max(getParam<Real>("pressure_max")),
175 652 : _num_T(getParam<unsigned int>("num_T")),
176 652 : _num_p(getParam<unsigned int>("num_p")),
177 826 : _fp(isParamValid("fp") ? &getUserObject<SinglePhaseFluidProperties>("fp") : nullptr),
178 652 : _allow_fp_and_tabulation(getParam<bool>("allow_fp_and_tabulation")),
179 978 : _interpolated_properties_enum(getParam<MultiMooseEnum>("interpolated_properties")),
180 : _interpolated_properties(),
181 326 : _interpolate_density(false),
182 326 : _interpolate_enthalpy(false),
183 326 : _interpolate_internal_energy(false),
184 326 : _interpolate_viscosity(false),
185 326 : _interpolate_k(false),
186 326 : _interpolate_c(false),
187 326 : _interpolate_cp(false),
188 326 : _interpolate_cv(false),
189 326 : _interpolate_entropy(false),
190 326 : _interpolate_pressure(false),
191 326 : _interpolate_temperature(false),
192 326 : _density_idx(libMesh::invalid_uint),
193 326 : _enthalpy_idx(libMesh::invalid_uint),
194 326 : _internal_energy_idx(libMesh::invalid_uint),
195 326 : _viscosity_idx(libMesh::invalid_uint),
196 326 : _k_idx(libMesh::invalid_uint),
197 326 : _c_idx(libMesh::invalid_uint),
198 326 : _cp_idx(libMesh::invalid_uint),
199 326 : _cv_idx(libMesh::invalid_uint),
200 326 : _entropy_idx(libMesh::invalid_uint),
201 326 : _p_idx(libMesh::invalid_uint),
202 326 : _T_idx(libMesh::invalid_uint),
203 326 : _csv_reader(_file_name_in, &_communicator),
204 652 : _construct_pT_from_ve(getParam<bool>("construct_pT_from_ve")),
205 652 : _construct_pT_from_vh(getParam<bool>("construct_pT_from_vh")),
206 326 : _initial_setup_done(false),
207 652 : _num_v(getParam<unsigned int>("num_v")),
208 652 : _num_e(getParam<unsigned int>("num_e")),
209 652 : _log_space_v(getParam<bool>("use_log_grid_v")),
210 652 : _log_space_e(getParam<bool>("use_log_grid_e")),
211 652 : _log_space_h(getParam<bool>("use_log_grid_h")),
212 1304 : _OOBBehavior(getParam<MooseEnum>("out_of_bounds_behavior"))
213 : {
214 : // Check that initial guess (used in Newton Method) is within min and max values
215 326 : checkInitialGuess();
216 : // Sanity check on minimum and maximum temperatures and pressures
217 322 : if (_temperature_max <= _temperature_min)
218 2 : mooseError("temperature_max must be greater than temperature_min");
219 320 : if (_pressure_max <= _pressure_min)
220 2 : mooseError("pressure_max must be greater than pressure_min");
221 :
222 : // Set (v,e) bounds if specified by the user
223 640 : if (isParamValid("e_min") && isParamValid("e_max"))
224 : {
225 0 : _e_min = getParam<Real>("e_min");
226 0 : _e_max = getParam<Real>("e_max");
227 0 : _e_bounds_specified = true;
228 : }
229 1268 : else if (isParamValid("e_min") || isParamValid("e_max"))
230 2 : paramError("e_min",
231 : "Either both or none of the min and max values of the specific internal energy "
232 : "should be specified");
233 : else
234 316 : _e_bounds_specified = false;
235 640 : if (isParamValid("v_min") && isParamValid("v_max"))
236 : {
237 4 : _v_min = getParam<Real>("v_min");
238 4 : _v_max = getParam<Real>("v_max");
239 2 : _v_bounds_specified = true;
240 : }
241 1252 : else if (isParamValid("v_min") || isParamValid("v_max"))
242 2 : paramError("v_min",
243 : "Either both or none of the min and max values of the specific volume "
244 : "should be specified");
245 : else
246 312 : _v_bounds_specified = false;
247 :
248 : // Handle out of bounds behavior parameters and deprecation
249 646 : if (isParamValid("error_on_out_of_bounds") && getParam<bool>("error_on_out_of_bounds") &&
250 6 : _OOBBehavior != Throw)
251 0 : paramError("out_of_bounds_behavior", "Inconsistent selection of out of bounds behavior.");
252 646 : else if (isParamValid("error_on_out_of_bounds") && !getParam<bool>("error_on_out_of_bounds"))
253 0 : _OOBBehavior = SetToClosestBound;
254 :
255 : // Lines starting with # in the data file are treated as comments
256 314 : _csv_reader.setComment("#");
257 :
258 : // Can only and must receive one source of data between fp and tabulations
259 314 : if (_fp && (!_file_name_in.empty() || !_file_name_ve_in.empty()) && !_allow_fp_and_tabulation)
260 0 : paramError("fluid_property_file",
261 : "Cannot supply both a fluid properties object with 'fp' and a source tabulation "
262 : "file with 'fluid_property_file', unless 'allow_fp_and_tabulation' is set to true");
263 314 : if (!_fp && _file_name_in.empty() && _file_name_ve_in.empty())
264 0 : paramError("fluid_property_file",
265 : "Either a fluid properties object with the parameter 'fp' and a source tabulation "
266 : "file with the parameter 'fluid_property_file' or 'fluid_property_ve_file' should "
267 : "be provided.");
268 314 : if (!_create_direct_pT_interpolations && !_create_direct_ve_interpolations)
269 0 : paramError("create_pT_interpolations", "Must create either (p,T) or (v,e) interpolations");
270 :
271 : // Some parameters are not used when reading a tabulation
272 434 : if (!_fp && !_file_name_in.empty() &&
273 806 : (isParamSetByUser("pressure_min") || isParamSetByUser("pressure_max") ||
274 650 : isParamSetByUser("temperature_min") || isParamSetByUser("temperature_max")))
275 36 : mooseWarning("User-specified bounds in pressure and temperature are ignored when reading a "
276 : "'fluid_property_file'. The tabulation bounds are selected "
277 : "from the bounds of the input tabulation.");
278 794 : if (!_fp && !_file_name_in.empty() && (isParamSetByUser("num_p") || isParamSetByUser("num_T")))
279 0 : mooseWarning("User-specified grid sizes in pressure and temperature are ignored when reading a "
280 : "'fluid_property_file'. The tabulation bounds are selected "
281 : "from the bounds of the input tabulation.");
282 342 : if (!_fp && !_file_name_ve_in.empty() &&
283 474 : (isParamSetByUser("v_min") || isParamSetByUser("v_max") || isParamSetByUser("e_min") ||
284 366 : isParamSetByUser("e_max")))
285 2 : mooseWarning(
286 : "User-specified bounds in specific volume and internal energy are ignored when reading a "
287 : "'fluid_property_ve_file'. The tabulation bounds are selected "
288 : "from the bounds of the input tabulation.");
289 416 : if (!_fp && !_file_name_ve_in.empty() && (isParamSetByUser("num_e") || isParamSetByUser("num_v")))
290 2 : mooseWarning("User-specified grid sizes in specific volume and internal energy are ignored "
291 : "when reading a 'fluid_property_ve_file'. The tabulation widths are read "
292 : "from the input tabulation.");
293 310 : if (!_file_name_ve_in.empty() && (_log_space_v || _log_space_e))
294 2 : mooseWarning(
295 : "User specfied logarithmic grids in specific volume and energy are ignored when reading a "
296 : "'fluid_properties_ve_file'. The tabulation grid is read from the input tabulation");
297 308 : }
298 :
299 : void
300 193 : TabulatedFluidProperties::initialSetup()
301 : {
302 193 : if (_initial_setup_done)
303 : return;
304 193 : _initial_setup_done = true;
305 :
306 193 : if (_create_direct_pT_interpolations)
307 : {
308 : // If the user specified a (p, T) tabulation to read, use that
309 169 : if (!_file_name_in.empty())
310 43 : readFileTabulationData(true);
311 : else
312 : {
313 126 : if (!_fp)
314 0 : paramError(
315 : "create_pT_interpolations",
316 : "No FluidProperties (specified with 'fp' parameter) exists. Either specify a 'fp' or "
317 : "specify a (p, T) tabulation file with the 'fluid_property_file' parameter");
318 126 : _console << name() + ": Generating (p, T) tabulated data\n";
319 126 : _console << std::flush;
320 :
321 126 : generateTabulatedData();
322 : }
323 : }
324 :
325 188 : if (_create_direct_ve_interpolations)
326 : {
327 : // If the user specified a (v, e) tabulation to read, use that
328 24 : if (!_file_name_ve_in.empty())
329 22 : readFileTabulationData(false);
330 : else
331 : {
332 2 : if (!_fp)
333 0 : paramError(
334 : "create_ve_interpolations",
335 : "No FluidProperties (specified with 'fp' parameter) exists. Either specify a 'fp' or "
336 : "specify a (v, e) tabulation file with the 'fluid_property_ve_file' parameter");
337 2 : _console << name() + ": Generating (v, e) tabulated data\n";
338 2 : _console << std::flush;
339 :
340 2 : generateVETabulatedData();
341 : }
342 : }
343 :
344 188 : computePropertyIndicesInInterpolationVectors();
345 188 : constructInterpolation();
346 :
347 : // Write tabulated data to file
348 188 : if (_save_file)
349 : {
350 24 : _console << name() + ": Writing tabulated data to " << _file_name_out << "\n";
351 48 : writeTabulatedData(_file_name_out);
352 : }
353 : }
354 :
355 : std::string
356 23 : TabulatedFluidProperties::fluidName() const
357 : {
358 23 : if (_fp)
359 23 : return _fp->fluidName();
360 : else
361 0 : return "TabulationFromFile";
362 : }
363 :
364 : Real
365 23 : TabulatedFluidProperties::molarMass() const
366 : {
367 23 : if (_fp)
368 23 : return _fp->molarMass();
369 : else
370 0 : FluidPropertiesForwardError("molarMass");
371 : }
372 :
373 : Real
374 68 : TabulatedFluidProperties::v_from_p_T(Real pressure, Real temperature) const
375 : {
376 68 : if (_interpolate_density && _create_direct_pT_interpolations)
377 : {
378 60 : checkInputVariables(pressure, temperature);
379 60 : return 1.0 / _property_ipol[_density_idx]->sample(pressure, temperature);
380 : }
381 : else
382 : {
383 8 : if (_fp)
384 8 : return 1.0 / _fp->rho_from_p_T(pressure, temperature);
385 : else
386 0 : paramError("fp", "No fluid properties or csv data provided for density.");
387 : }
388 : }
389 :
390 : void
391 18495375 : TabulatedFluidProperties::v_from_p_T(
392 : Real pressure, Real temperature, Real & v, Real & dv_dp, Real & dv_dT) const
393 : {
394 18495375 : Real rho = 0, drho_dp = 0, drho_dT = 0;
395 18495375 : if (_interpolate_density && _create_direct_pT_interpolations)
396 : {
397 18495375 : checkInputVariables(pressure, temperature);
398 18495375 : _property_ipol[_density_idx]->sampleValueAndDerivatives(
399 : pressure, temperature, rho, drho_dp, drho_dT);
400 : }
401 : else
402 : {
403 0 : if (_fp)
404 0 : _fp->rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
405 : else
406 0 : paramError("fp", "No fluid properties or csv data provided for density.");
407 : }
408 : // convert from rho to v
409 18495375 : v = 1.0 / rho;
410 18495375 : dv_dp = -drho_dp / (rho * rho);
411 18495375 : dv_dT = -drho_dT / (rho * rho);
412 18495375 : }
413 :
414 : Real
415 702 : TabulatedFluidProperties::rho_from_p_T(Real pressure, Real temperature) const
416 : {
417 702 : if (_interpolate_density && _create_direct_pT_interpolations)
418 : {
419 699 : checkInputVariables(pressure, temperature);
420 695 : return _property_ipol[_density_idx]->sample(pressure, temperature);
421 : }
422 : else
423 : {
424 3 : if (_fp)
425 3 : return _fp->rho_from_p_T(pressure, temperature);
426 : else
427 0 : paramError("fp", "No fluid properties or csv data provided for density.");
428 : }
429 : }
430 :
431 : void
432 55 : TabulatedFluidProperties::rho_from_p_T(
433 : Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
434 : {
435 55 : if (_interpolate_density && _create_direct_pT_interpolations)
436 : {
437 52 : checkInputVariables(pressure, temperature);
438 52 : _property_ipol[_density_idx]->sampleValueAndDerivatives(
439 : pressure, temperature, rho, drho_dp, drho_dT);
440 : }
441 : else
442 : {
443 3 : if (_fp)
444 3 : _fp->rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
445 : else
446 0 : paramError("fp", "No fluid properties or csv data provided for density.");
447 : }
448 55 : }
449 :
450 : void
451 0 : TabulatedFluidProperties::rho_from_p_T(const ADReal & pressure,
452 : const ADReal & temperature,
453 : ADReal & rho,
454 : ADReal & drho_dp,
455 : ADReal & drho_dT) const
456 : {
457 0 : if (_interpolate_density && _create_direct_pT_interpolations)
458 : {
459 0 : ADReal p = pressure, T = temperature;
460 0 : checkInputVariables(p, T);
461 0 : _property_ipol[_density_idx]->sampleValueAndDerivatives(p, T, rho, drho_dp, drho_dT);
462 0 : }
463 : else
464 : {
465 0 : if (_fp)
466 0 : _fp->rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
467 : else
468 0 : paramError("fp", "No fluid properties or csv data provided for density.");
469 : }
470 0 : }
471 :
472 : Real
473 6 : TabulatedFluidProperties::rho_from_p_s(Real p, Real s) const
474 : {
475 6 : Real T = T_from_p_s(p, s);
476 6 : return rho_from_p_T(p, T);
477 : }
478 :
479 : void
480 1 : TabulatedFluidProperties::rho_from_p_s(
481 : Real p, Real s, Real & rho, Real & drho_dp, Real & drho_ds) const
482 : {
483 : Real T, dT_dp, dT_ds;
484 1 : T_from_p_s(p, s, T, dT_dp, dT_ds);
485 : Real drho_dp_T, drho_dT;
486 1 : rho_from_p_T(p, T, rho, drho_dp_T, drho_dT);
487 1 : drho_dp = drho_dT * dT_dp + drho_dp_T;
488 1 : drho_ds = drho_dT * dT_ds;
489 1 : }
490 :
491 : Real
492 330 : TabulatedFluidProperties::e_from_p_T(Real pressure, Real temperature) const
493 : {
494 330 : if (_interpolate_internal_energy && _create_direct_pT_interpolations)
495 : {
496 320 : checkInputVariables(pressure, temperature);
497 320 : return _property_ipol[_internal_energy_idx]->sample(pressure, temperature);
498 : }
499 : else
500 : {
501 10 : if (_fp)
502 10 : return _fp->e_from_p_T(pressure, temperature);
503 : else
504 0 : paramError("fp", "No fluid properties or csv data provided for internal energy.");
505 : }
506 : }
507 :
508 : void
509 9142340 : TabulatedFluidProperties::e_from_p_T(
510 : Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const
511 : {
512 9142340 : if (_interpolate_internal_energy && _create_direct_pT_interpolations)
513 : {
514 9142338 : checkInputVariables(pressure, temperature);
515 9142338 : _property_ipol[_internal_energy_idx]->sampleValueAndDerivatives(
516 : pressure, temperature, e, de_dp, de_dT);
517 : }
518 : else
519 : {
520 2 : if (_fp)
521 2 : _fp->e_from_p_T(pressure, temperature, e, de_dp, de_dT);
522 : else
523 0 : paramError("fp", "No fluid properties or csv data provided for internal energy.");
524 : }
525 9142340 : }
526 :
527 : Real
528 5 : TabulatedFluidProperties::e_from_p_rho(Real pressure, Real rho) const
529 : {
530 5 : Real T = T_from_p_rho(pressure, rho);
531 5 : Real e = e_from_p_T(pressure, T);
532 5 : return e;
533 : }
534 :
535 : void
536 1 : TabulatedFluidProperties::e_from_p_rho(
537 : Real pressure, Real rho, Real & e, Real & de_dp, Real & de_drho) const
538 : {
539 : // get derivatives of T wrt to pressure and density
540 : Real T, dT_dp, dT_drho;
541 1 : T_from_p_rho(pressure, rho, T, dT_dp, dT_drho);
542 :
543 : // Get e, then derivatives of e wrt pressure and temperature
544 : Real de_dp_at_const_T, de_dT;
545 1 : e_from_p_T(pressure, T, e, de_dp_at_const_T, de_dT);
546 :
547 : // Get the derivatives of density wrt pressure and temperature
548 : Real rho_pT, drho_dp, drho_dT;
549 1 : rho_from_p_T(pressure, T, rho_pT, drho_dp, drho_dT);
550 :
551 : // derivatives of e wrt pressure and rho (what we want from e_from_p_rho)
552 1 : de_drho = de_dT * dT_drho;
553 1 : de_dp = de_dp_at_const_T - (de_drho * drho_dp);
554 1 : }
555 :
556 : Real
557 16 : TabulatedFluidProperties::T_from_p_rho(Real pressure, Real rho) const
558 : {
559 : auto lambda = [&](Real p, Real current_T, Real & new_rho, Real & drho_dp, Real & drho_dT)
560 48 : { rho_from_p_T(p, current_T, new_rho, drho_dp, drho_dT); };
561 16 : Real T = FluidPropertiesUtils::NewtonSolve(pressure,
562 : rho,
563 16 : _T_initial_guess,
564 16 : _tolerance,
565 : lambda,
566 16 : name() + "::T_from_p_rho",
567 16 : _max_newton_its)
568 16 : .first;
569 : // check for nans
570 16 : if (std::isnan(T))
571 0 : mooseError("Conversion from pressure (p = ",
572 : pressure,
573 : ") and density (rho = ",
574 : rho,
575 : ") to temperature failed to converge.");
576 16 : return T;
577 : }
578 :
579 : void
580 2 : TabulatedFluidProperties::T_from_p_rho(
581 : Real pressure, Real rho, Real & T, Real & dT_dp, Real & dT_drho) const
582 : {
583 2 : T = T_from_p_rho(pressure, rho);
584 : Real eps = 1e-8;
585 2 : dT_dp = (T_from_p_rho(pressure * (1 + eps), rho) - T) / (eps * pressure);
586 2 : dT_drho = (T_from_p_rho(pressure, rho * (1 + eps)) - T) / (eps * rho);
587 2 : }
588 :
589 : Real
590 18 : TabulatedFluidProperties::T_from_p_s(Real pressure, Real s) const
591 : {
592 : auto lambda = [&](Real p, Real current_T, Real & new_s, Real & ds_dp, Real & ds_dT)
593 50 : { s_from_p_T(p, current_T, new_s, ds_dp, ds_dT); };
594 18 : Real T = FluidPropertiesUtils::NewtonSolve(pressure,
595 : s,
596 18 : _T_initial_guess,
597 18 : _tolerance,
598 : lambda,
599 18 : name() + "::T_from_p_s",
600 18 : _max_newton_its)
601 18 : .first;
602 : // check for nans
603 18 : if (std::isnan(T))
604 0 : mooseError("Conversion from pressure (p = ",
605 : pressure,
606 : ") and entropy (s = ",
607 : s,
608 : ") to temperature failed to converge.");
609 18 : return T;
610 : }
611 :
612 : void
613 2 : TabulatedFluidProperties::T_from_p_s(
614 : Real pressure, Real s, Real & T, Real & dT_dp, Real & dT_ds) const
615 : {
616 2 : T = T_from_p_s(pressure, s);
617 : Real eps = 1e-8;
618 2 : dT_dp = (T_from_p_s(pressure * (1 + eps), s) - T) / (eps * pressure);
619 2 : dT_ds = (T_from_p_s(pressure, s * (1 + eps)) - T) / (eps * s);
620 2 : }
621 :
622 : Real
623 1563 : TabulatedFluidProperties::h_from_p_T(Real pressure, Real temperature) const
624 : {
625 1563 : if (_interpolate_enthalpy && _create_direct_pT_interpolations)
626 : {
627 1562 : checkInputVariables(pressure, temperature);
628 1562 : return _property_ipol[_enthalpy_idx]->sample(pressure, temperature);
629 : }
630 : else
631 : {
632 1 : if (_fp)
633 1 : return _fp->h_from_p_T(pressure, temperature);
634 : else
635 0 : paramError("fp", "No fluid properties or csv data provided for enthalpy.");
636 : }
637 : }
638 :
639 : ADReal
640 0 : TabulatedFluidProperties::h_from_p_T(const ADReal & pressure, const ADReal & temperature) const
641 : {
642 0 : if (_fp) // Assuming _fp can handle ADReal types
643 0 : return _fp->h_from_p_T(pressure, temperature);
644 : else
645 0 : FluidPropertiesForwardError("h_from_p_T");
646 : }
647 :
648 : void
649 9353048 : TabulatedFluidProperties::h_from_p_T(
650 : Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const
651 : {
652 9353048 : if (_interpolate_enthalpy && _create_direct_pT_interpolations)
653 : {
654 9353048 : checkInputVariables(pressure, temperature);
655 9353048 : _property_ipol[_enthalpy_idx]->sampleValueAndDerivatives(
656 : pressure, temperature, h, dh_dp, dh_dT);
657 : }
658 : else
659 : {
660 0 : if (_fp)
661 0 : _fp->h_from_p_T(pressure, temperature, h, dh_dp, dh_dT);
662 : else
663 0 : paramError("fp", "No fluid properties or csv data provided for enthalpy.");
664 : }
665 9353048 : }
666 :
667 : Real
668 1092 : TabulatedFluidProperties::mu_from_p_T(Real pressure, Real temperature) const
669 : {
670 1092 : if (_interpolate_viscosity && _create_direct_pT_interpolations)
671 : {
672 1090 : checkInputVariables(pressure, temperature);
673 1090 : return _property_ipol[_viscosity_idx]->sample(pressure, temperature);
674 : }
675 : else
676 : {
677 2 : if (_fp)
678 2 : return _fp->mu_from_p_T(pressure, temperature);
679 : else
680 0 : paramError("fp", "No fluid properties or csv data provided for viscosity.");
681 : }
682 : }
683 :
684 : void
685 7 : TabulatedFluidProperties::mu_from_p_T(
686 : Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
687 : {
688 7 : if (_interpolate_viscosity && _create_direct_pT_interpolations)
689 : {
690 5 : checkInputVariables(pressure, temperature);
691 5 : _property_ipol[_viscosity_idx]->sampleValueAndDerivatives(
692 : pressure, temperature, mu, dmu_dp, dmu_dT);
693 : }
694 : else
695 : {
696 2 : if (_fp)
697 2 : _fp->mu_from_p_T(pressure, temperature, mu, dmu_dp, dmu_dT);
698 : else
699 0 : paramError("fp", "No fluid properties or csv data provided for viscosity.");
700 : }
701 7 : }
702 :
703 : Real
704 1072 : TabulatedFluidProperties::c_from_p_T(Real pressure, Real temperature) const
705 : {
706 1072 : if (_interpolate_c && _create_direct_pT_interpolations)
707 : {
708 1072 : checkInputVariables(pressure, temperature);
709 1072 : return _property_ipol[_c_idx]->sample(pressure, temperature);
710 : }
711 : else
712 : {
713 0 : if (_fp)
714 0 : return _fp->c_from_p_T(pressure, temperature);
715 : else
716 0 : paramError("interpolated_properties", "No data to interpolate for speed of sound.");
717 : }
718 : }
719 :
720 : void
721 0 : TabulatedFluidProperties::c_from_p_T(
722 : Real pressure, Real temperature, Real & c, Real & dc_dp, Real & dc_dT) const
723 : {
724 0 : if (_interpolate_c && _create_direct_pT_interpolations)
725 : {
726 0 : checkInputVariables(pressure, temperature);
727 0 : _property_ipol[_c_idx]->sampleValueAndDerivatives(pressure, temperature, c, dc_dp, dc_dT);
728 : }
729 : else
730 : {
731 0 : if (_fp)
732 0 : _fp->c_from_p_T(pressure, temperature, c, dc_dp, dc_dT);
733 : else
734 0 : paramError("interpolated_properties", "No data to interpolate for speed of sound.");
735 : }
736 0 : }
737 :
738 : Real
739 1307 : TabulatedFluidProperties::cp_from_p_T(Real pressure, Real temperature) const
740 : {
741 1307 : if (_interpolate_cp && _create_direct_pT_interpolations)
742 : {
743 1089 : checkInputVariables(pressure, temperature);
744 1089 : return _property_ipol[_cp_idx]->sample(pressure, temperature);
745 : }
746 : else
747 : {
748 218 : if (_fp)
749 218 : return _fp->cp_from_p_T(pressure, temperature);
750 : else
751 0 : paramError("interpolated_properties",
752 : "No data to interpolate for specific heat capacity at constant pressure.");
753 : }
754 : }
755 :
756 : void
757 4 : TabulatedFluidProperties::cp_from_p_T(
758 : Real pressure, Real temperature, Real & cp, Real & dcp_dp, Real & dcp_dT) const
759 : {
760 4 : if (_interpolate_cp && _create_direct_pT_interpolations)
761 : {
762 4 : checkInputVariables(pressure, temperature);
763 4 : _property_ipol[_cp_idx]->sampleValueAndDerivatives(pressure, temperature, cp, dcp_dp, dcp_dT);
764 : }
765 : else
766 : {
767 0 : if (_fp)
768 0 : _fp->cp_from_p_T(pressure, temperature, cp, dcp_dp, dcp_dT);
769 : else
770 0 : paramError("interpolated_properties",
771 : "No data to interpolate for specific heat capacity at constant pressure.");
772 : }
773 4 : }
774 :
775 : Real
776 1091 : TabulatedFluidProperties::cv_from_p_T(Real pressure, Real temperature) const
777 : {
778 1091 : if (_interpolate_cv && _create_direct_pT_interpolations)
779 : {
780 1089 : checkInputVariables(pressure, temperature);
781 1089 : return _property_ipol[_cv_idx]->sample(pressure, temperature);
782 : }
783 : else
784 : {
785 2 : if (_fp)
786 2 : return _fp->cv_from_p_T(pressure, temperature);
787 : else
788 0 : paramError("interpolated_properties",
789 : "No data to interpolate for specific heat capacity at constant volume.");
790 : }
791 : }
792 :
793 : void
794 4 : TabulatedFluidProperties::cv_from_p_T(
795 : Real pressure, Real temperature, Real & cv, Real & dcv_dp, Real & dcv_dT) const
796 : {
797 4 : if (_interpolate_cv)
798 : {
799 4 : checkInputVariables(pressure, temperature);
800 4 : _property_ipol[_cv_idx]->sampleValueAndDerivatives(pressure, temperature, cv, dcv_dp, dcv_dT);
801 : }
802 : else
803 : {
804 0 : if (_fp)
805 0 : _fp->cv_from_p_T(pressure, temperature, cv, dcv_dp, dcv_dT);
806 : else
807 0 : paramError("interpolated_properties",
808 : "No data to interpolate for specific heat capacity at constant volume.");
809 : }
810 4 : }
811 :
812 : Real
813 1091 : TabulatedFluidProperties::k_from_p_T(Real pressure, Real temperature) const
814 : {
815 1091 : if (_interpolate_k && _create_direct_pT_interpolations)
816 : {
817 1089 : checkInputVariables(pressure, temperature);
818 1089 : return _property_ipol[_k_idx]->sample(pressure, temperature);
819 : }
820 : else
821 : {
822 2 : if (_fp)
823 2 : return _fp->k_from_p_T(pressure, temperature);
824 : else
825 0 : paramError("interpolated_properties", "No data to interpolate for thermal conductivity.");
826 : }
827 : }
828 :
829 : void
830 4 : TabulatedFluidProperties::k_from_p_T(
831 : Real pressure, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
832 : {
833 4 : if (_interpolate_k && _create_direct_pT_interpolations)
834 : {
835 4 : checkInputVariables(pressure, temperature);
836 4 : return _property_ipol[_k_idx]->sampleValueAndDerivatives(
837 4 : pressure, temperature, k, dk_dp, dk_dT);
838 : }
839 : else
840 : {
841 0 : if (_fp)
842 0 : return _fp->k_from_p_T(pressure, temperature, k, dk_dp, dk_dT);
843 : else
844 0 : paramError("interpolated_properties", "No data to interpolate for thermal conductivity.");
845 : }
846 : }
847 :
848 : Real
849 1101 : TabulatedFluidProperties::s_from_p_T(Real pressure, Real temperature) const
850 : {
851 1101 : if (_interpolate_entropy && _create_direct_pT_interpolations)
852 : {
853 1099 : checkInputVariables(pressure, temperature);
854 1099 : return _property_ipol[_entropy_idx]->sample(pressure, temperature);
855 : }
856 : else
857 : {
858 2 : if (_fp)
859 2 : return _fp->s_from_p_T(pressure, temperature);
860 : else
861 0 : paramError("interpolated_properties", "No data to interpolate for entropy.");
862 : }
863 : }
864 :
865 : void
866 50 : TabulatedFluidProperties::s_from_p_T(Real p, Real T, Real & s, Real & ds_dp, Real & ds_dT) const
867 : {
868 50 : if (_interpolate_entropy && _create_direct_pT_interpolations)
869 : {
870 50 : checkInputVariables(p, T);
871 50 : _property_ipol[_entropy_idx]->sampleValueAndDerivatives(p, T, s, ds_dp, ds_dT);
872 : }
873 : else
874 : {
875 0 : if (_fp)
876 0 : _fp->s_from_p_T(p, T, s, ds_dp, ds_dT);
877 : else
878 0 : paramError("interpolated_properties", "No data to interpolate for entropy.");
879 : }
880 50 : }
881 :
882 : Real
883 15 : TabulatedFluidProperties::e_from_v_h(Real v, Real h) const
884 : {
885 15 : if (_construct_pT_from_vh)
886 : {
887 4 : const Real p = _p_from_v_h_ipol->sample(v, h);
888 4 : const Real T = _T_from_v_h_ipol->sample(v, h);
889 4 : return e_from_p_T(p, T);
890 : }
891 11 : else if (_create_direct_ve_interpolations)
892 : {
893 : // Lambda computes h from v and the current_e
894 : auto lambda = [&](Real v, Real current_e, Real & new_h, Real & dh_dv, Real & dh_de)
895 33 : { h_from_v_e(v, current_e, new_h, dh_dv, dh_de); };
896 11 : Real e = FluidPropertiesUtils::NewtonSolve(v,
897 : h,
898 11 : /*e initial guess*/ h - _p_initial_guess * v,
899 11 : _tolerance,
900 : lambda,
901 22 : name() + "::e_from_v_h",
902 11 : _max_newton_its)
903 11 : .first;
904 : return e;
905 : }
906 0 : else if (_fp)
907 0 : return _fp->e_from_v_h(v, h);
908 : else
909 0 : mooseError(__PRETTY_FUNCTION__,
910 : "\nNo tabulation or fluid property 'fp' object to compute value");
911 : }
912 :
913 : void
914 3 : TabulatedFluidProperties::e_from_v_h(Real v, Real h, Real & e, Real & de_dv, Real & de_dh) const
915 : {
916 3 : if (_construct_pT_from_vh)
917 : {
918 1 : Real p = 0, dp_dv = 0, dp_dh = 0;
919 1 : _p_from_v_h_ipol->sampleValueAndDerivatives(v, h, p, dp_dv, dp_dh);
920 1 : Real T = 0, dT_dv = 0, dT_dh = 0;
921 1 : _T_from_v_h_ipol->sampleValueAndDerivatives(v, h, T, dT_dv, dT_dh);
922 : Real de_dp, de_dT;
923 1 : e_from_p_T(p, T, e, de_dp, de_dT);
924 1 : de_dv = de_dp * dp_dv + de_dT * dT_dv;
925 1 : de_dh = de_dp * dp_dh + de_dT * dT_dh;
926 : }
927 2 : else if (_create_direct_ve_interpolations)
928 : {
929 : // Lambda computes h from v and the current_e
930 : auto lambda = [&](Real v, Real current_e, Real & new_h, Real & dh_dv, Real & dh_de)
931 6 : { h_from_v_e(v, current_e, new_h, dh_dv, dh_de); };
932 : const auto e_data =
933 2 : FluidPropertiesUtils::NewtonSolve(v,
934 : h,
935 2 : /*e initial guess*/ h - _p_initial_guess * v,
936 2 : _tolerance,
937 : lambda,
938 2 : name() + "::e_from_v_h",
939 2 : _max_newton_its);
940 2 : e = e_data.first;
941 : // Finite difference approximation
942 2 : const auto e2 = e_from_v_h(v * (1 + TOLERANCE), h);
943 2 : de_dv = (e2 - e) / (TOLERANCE * v);
944 2 : de_dh = 1. / e_data.second;
945 : }
946 0 : else if (_fp)
947 0 : _fp->e_from_v_h(v, h, e, de_dv, de_dh);
948 : else
949 0 : mooseError(__PRETTY_FUNCTION__,
950 : "\nNo tabulation or fluid property 'fp' object to compute value");
951 3 : }
952 :
953 : std::vector<Real>
954 3 : TabulatedFluidProperties::henryCoefficients() const
955 : {
956 3 : if (_fp)
957 3 : return _fp->henryCoefficients();
958 : else
959 0 : FluidPropertiesForwardError("henryCoefficients");
960 : }
961 :
962 : Real
963 1 : TabulatedFluidProperties::vaporPressure(Real temperature) const
964 : {
965 1 : if (_fp)
966 1 : return _fp->vaporPressure(temperature);
967 : else
968 0 : FluidPropertiesForwardError("vaporPressure");
969 : }
970 :
971 : void
972 0 : TabulatedFluidProperties::vaporPressure(Real temperature, Real & psat, Real & dpsat_dT) const
973 : {
974 0 : if (_fp)
975 0 : _fp->vaporPressure(temperature, psat, dpsat_dT);
976 : else
977 0 : FluidPropertiesForwardError("vaporPressure");
978 0 : }
979 :
980 : Real
981 0 : TabulatedFluidProperties::vaporTemperature(Real pressure) const
982 : {
983 0 : if (_fp)
984 0 : return _fp->vaporTemperature(pressure);
985 : else
986 0 : FluidPropertiesForwardError("vaporTemperature");
987 : }
988 :
989 : void
990 0 : TabulatedFluidProperties::vaporTemperature(Real pressure, Real & Tsat, Real & dTsat_dp) const
991 : {
992 :
993 0 : if (_fp)
994 0 : _fp->vaporTemperature(pressure, Tsat, dTsat_dp);
995 : else
996 0 : FluidPropertiesForwardError("vaporTemperature");
997 0 : }
998 :
999 : Real
1000 1 : TabulatedFluidProperties::triplePointPressure() const
1001 : {
1002 :
1003 1 : if (_fp)
1004 1 : return _fp->triplePointPressure();
1005 : else
1006 0 : FluidPropertiesForwardError("triplePointPressure");
1007 : }
1008 :
1009 : Real
1010 1 : TabulatedFluidProperties::triplePointTemperature() const
1011 : {
1012 :
1013 1 : if (_fp)
1014 1 : return _fp->triplePointTemperature();
1015 : else
1016 0 : FluidPropertiesForwardError("triplePointTemperature");
1017 : }
1018 :
1019 : Real
1020 1 : TabulatedFluidProperties::criticalPressure() const
1021 : {
1022 :
1023 1 : if (_fp)
1024 1 : return _fp->criticalPressure();
1025 : else
1026 0 : FluidPropertiesForwardError("criticalPressure");
1027 : }
1028 :
1029 : Real
1030 1 : TabulatedFluidProperties::criticalTemperature() const
1031 : {
1032 :
1033 1 : if (_fp)
1034 1 : return _fp->criticalTemperature();
1035 : else
1036 0 : FluidPropertiesForwardError("criticalTemperature");
1037 : }
1038 :
1039 : Real
1040 1 : TabulatedFluidProperties::criticalDensity() const
1041 : {
1042 1 : if (_fp)
1043 1 : return _fp->criticalDensity();
1044 : else
1045 0 : FluidPropertiesForwardError("criticalDensity");
1046 : }
1047 :
1048 : Real
1049 46 : TabulatedFluidProperties::p_from_v_e(Real v, Real e) const
1050 : {
1051 46 : if (_interpolate_pressure && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1052 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1053 46 : checkInputVariablesVE(v, e);
1054 :
1055 46 : if (_create_direct_ve_interpolations && _interpolate_pressure)
1056 15 : return _property_ve_ipol[_p_idx]->sample(v, e);
1057 31 : else if (_construct_pT_from_ve)
1058 31 : return _p_from_v_e_ipol->sample(v, e);
1059 0 : else if (_fp)
1060 0 : return _fp->p_from_v_e(v, e);
1061 : else
1062 0 : mooseError(__PRETTY_FUNCTION__,
1063 : "\nNo tabulation or fluid property 'fp' object to compute value");
1064 : }
1065 :
1066 : void
1067 47 : TabulatedFluidProperties::p_from_v_e(Real v, Real e, Real & p, Real & dp_dv, Real & dp_de) const
1068 : {
1069 47 : if (_interpolate_pressure && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1070 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1071 47 : checkInputVariablesVE(v, e);
1072 :
1073 47 : if (_create_direct_ve_interpolations && _interpolate_pressure)
1074 43 : _property_ve_ipol[_p_idx]->sampleValueAndDerivatives(v, e, p, dp_dv, dp_de);
1075 4 : else if (_construct_pT_from_ve)
1076 4 : _p_from_v_e_ipol->sampleValueAndDerivatives(v, e, p, dp_dv, dp_de);
1077 0 : else if (_fp)
1078 0 : _fp->p_from_v_e(v, e, p, dp_dv, dp_de);
1079 : else
1080 0 : mooseError(__PRETTY_FUNCTION__,
1081 : "\nNo tabulation or fluid property 'fp' object to compute value");
1082 47 : }
1083 :
1084 : Real
1085 46 : TabulatedFluidProperties::T_from_v_e(Real v, Real e) const
1086 : {
1087 46 : if (_interpolate_temperature && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1088 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1089 46 : checkInputVariablesVE(v, e);
1090 :
1091 46 : if (_create_direct_ve_interpolations && _interpolate_temperature)
1092 15 : return _property_ve_ipol[_T_idx]->sample(v, e);
1093 31 : else if (_construct_pT_from_ve)
1094 31 : return _T_from_v_e_ipol->sample(v, e);
1095 0 : else if (_fp)
1096 0 : return _fp->T_from_v_e(v, e);
1097 : else
1098 0 : mooseError(__PRETTY_FUNCTION__,
1099 : "\nNo tabulation or fluid property 'fp' object to compute value");
1100 : }
1101 :
1102 : void
1103 8 : TabulatedFluidProperties::T_from_v_e(Real v, Real e, Real & T, Real & dT_dv, Real & dT_de) const
1104 : {
1105 8 : if (_interpolate_temperature && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1106 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1107 8 : checkInputVariablesVE(v, e);
1108 :
1109 8 : if (_create_direct_ve_interpolations && _interpolate_temperature)
1110 4 : _property_ve_ipol[_T_idx]->sampleValueAndDerivatives(v, e, T, dT_dv, dT_de);
1111 4 : else if (_construct_pT_from_ve)
1112 4 : _T_from_v_e_ipol->sampleValueAndDerivatives(v, e, T, dT_dv, dT_de);
1113 0 : else if (_fp)
1114 0 : _fp->T_from_v_e(v, e, T, dT_dv, dT_de);
1115 : else
1116 0 : mooseError(__PRETTY_FUNCTION__,
1117 : "\nNo tabulation or fluid property 'fp' object to compute value");
1118 8 : }
1119 :
1120 : Real
1121 835 : TabulatedFluidProperties::c_from_v_e(Real v, Real e) const
1122 : {
1123 835 : if (_interpolate_c && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1124 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1125 835 : checkInputVariablesVE(v, e);
1126 :
1127 835 : if (_create_direct_ve_interpolations && _interpolate_c)
1128 11 : return _property_ve_ipol[_c_idx]->sample(v, e);
1129 824 : else if (_construct_pT_from_ve)
1130 : {
1131 824 : Real p = _p_from_v_e_ipol->sample(v, e);
1132 824 : Real T = _T_from_v_e_ipol->sample(v, e);
1133 824 : return c_from_p_T(p, T);
1134 : }
1135 0 : else if (_fp)
1136 0 : return _fp->c_from_v_e(v, e);
1137 : else
1138 0 : mooseError(__PRETTY_FUNCTION__,
1139 : "\nNo tabulation or fluid property 'fp' object to compute value");
1140 : }
1141 :
1142 : void
1143 1 : TabulatedFluidProperties::c_from_v_e(Real v, Real e, Real & c, Real & dc_dv, Real & dc_de) const
1144 : {
1145 1 : if (_interpolate_c && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1146 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1147 1 : checkInputVariablesVE(v, e);
1148 :
1149 1 : if (_create_direct_ve_interpolations && _interpolate_c)
1150 1 : _property_ve_ipol[_c_idx]->sampleValueAndDerivatives(v, e, c, dc_dv, dc_de);
1151 0 : else if (_construct_pT_from_ve)
1152 : {
1153 : Real p, dp_dv, dp_de;
1154 0 : _p_from_v_e_ipol->sampleValueAndDerivatives(v, e, p, dp_dv, dp_de);
1155 : Real T, dT_dv, dT_de;
1156 0 : _T_from_v_e_ipol->sampleValueAndDerivatives(v, e, T, dT_dv, dT_de);
1157 : Real dc_dp, dc_dT;
1158 0 : c_from_p_T(p, T, c, dc_dp, dc_dT);
1159 0 : dc_dv = dc_dp * dp_dv + dc_dT * dT_dv;
1160 0 : dc_de = dc_dp * dp_de + dc_dT * dT_de;
1161 : }
1162 0 : else if (_fp)
1163 0 : _fp->c_from_v_e(v, e, c, dc_dv, dc_de);
1164 : else
1165 0 : mooseError(__PRETTY_FUNCTION__,
1166 : "\nNo tabulation or fluid property 'fp' object to compute value");
1167 1 : }
1168 :
1169 : Real
1170 853 : TabulatedFluidProperties::cp_from_v_e(Real v, Real e) const
1171 : {
1172 853 : if (_interpolate_cp && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1173 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1174 853 : checkInputVariablesVE(v, e);
1175 :
1176 853 : if (_create_direct_ve_interpolations && _interpolate_cp)
1177 19 : return _property_ve_ipol[_cp_idx]->sample(v, e);
1178 834 : else if (_construct_pT_from_ve)
1179 : {
1180 834 : Real p = _p_from_v_e_ipol->sample(v, e);
1181 834 : Real T = _T_from_v_e_ipol->sample(v, e);
1182 834 : return cp_from_p_T(p, T);
1183 : }
1184 0 : else if (_fp)
1185 0 : return _fp->cp_from_v_e(v, e);
1186 : else
1187 0 : mooseError(__PRETTY_FUNCTION__,
1188 : "\nNo tabulation or fluid property 'fp' object to compute value");
1189 : }
1190 :
1191 : void
1192 6 : TabulatedFluidProperties::cp_from_v_e(Real v, Real e, Real & cp, Real & dcp_dv, Real & dcp_de) const
1193 : {
1194 6 : if (_interpolate_cp && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1195 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1196 6 : checkInputVariablesVE(v, e);
1197 :
1198 6 : if (_create_direct_ve_interpolations && _interpolate_cp)
1199 3 : _property_ve_ipol[_cp_idx]->sampleValueAndDerivatives(v, e, cp, dcp_dv, dcp_de);
1200 3 : else if (_construct_pT_from_ve)
1201 : {
1202 : Real p, dp_dv, dp_de;
1203 3 : _p_from_v_e_ipol->sampleValueAndDerivatives(v, e, p, dp_dv, dp_de);
1204 : Real T, dT_dv, dT_de;
1205 3 : _T_from_v_e_ipol->sampleValueAndDerivatives(v, e, T, dT_dv, dT_de);
1206 : Real dcp_dp, dcp_dT;
1207 3 : cp_from_p_T(p, T, cp, dcp_dp, dcp_dT);
1208 3 : dcp_dv = dcp_dp * dp_dv + dcp_dT * dT_dv;
1209 3 : dcp_de = dcp_dp * dp_de + dcp_dT * dT_de;
1210 : }
1211 0 : else if (_fp)
1212 0 : _fp->cp_from_v_e(v, e, cp, dcp_dv, dcp_de);
1213 : else
1214 0 : mooseError(__PRETTY_FUNCTION__,
1215 : "\nNo tabulation or fluid property 'fp' object to compute value");
1216 6 : }
1217 :
1218 : Real
1219 853 : TabulatedFluidProperties::cv_from_v_e(Real v, Real e) const
1220 : {
1221 853 : if (_interpolate_cv && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1222 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1223 853 : checkInputVariablesVE(v, e);
1224 :
1225 853 : if (_create_direct_ve_interpolations && _interpolate_cv)
1226 19 : return _property_ve_ipol[_cv_idx]->sample(v, e);
1227 834 : else if (_construct_pT_from_ve)
1228 : {
1229 834 : Real p = _p_from_v_e_ipol->sample(v, e);
1230 834 : Real T = _T_from_v_e_ipol->sample(v, e);
1231 834 : return cv_from_p_T(p, T);
1232 : }
1233 0 : else if (_fp)
1234 0 : return _fp->cv_from_v_e(v, e);
1235 : else
1236 0 : mooseError(__PRETTY_FUNCTION__,
1237 : "\nNo tabulation or fluid property 'fp' object to compute value");
1238 : }
1239 :
1240 : void
1241 6 : TabulatedFluidProperties::cv_from_v_e(Real v, Real e, Real & cv, Real & dcv_dv, Real & dcv_de) const
1242 : {
1243 6 : if (_interpolate_cv && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1244 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1245 6 : checkInputVariablesVE(v, e);
1246 :
1247 6 : if (_create_direct_ve_interpolations && _interpolate_cv)
1248 3 : _property_ve_ipol[_cv_idx]->sampleValueAndDerivatives(v, e, cv, dcv_dv, dcv_de);
1249 3 : else if (_construct_pT_from_ve)
1250 : {
1251 : Real p, dp_dv, dp_de;
1252 3 : _p_from_v_e_ipol->sampleValueAndDerivatives(v, e, p, dp_dv, dp_de);
1253 : Real T, dT_dv, dT_de;
1254 3 : _T_from_v_e_ipol->sampleValueAndDerivatives(v, e, T, dT_dv, dT_de);
1255 : Real dcv_dp, dcv_dT;
1256 3 : cv_from_p_T(p, T, cv, dcv_dp, dcv_dT);
1257 3 : dcv_dv = dcv_dp * dp_dv + dcv_dT * dT_dv;
1258 3 : dcv_de = dcv_dp * dp_de + dcv_dT * dT_de;
1259 : }
1260 0 : else if (_fp)
1261 0 : _fp->cv_from_v_e(v, e, cv, dcv_dv, dcv_de);
1262 : else
1263 0 : mooseError(__PRETTY_FUNCTION__,
1264 : "\nNo tabulation or fluid property 'fp' object to compute value");
1265 6 : }
1266 :
1267 : Real
1268 853 : TabulatedFluidProperties::mu_from_v_e(Real v, Real e) const
1269 : {
1270 853 : if (_interpolate_viscosity && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1271 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1272 853 : checkInputVariablesVE(v, e);
1273 :
1274 853 : if (_create_direct_ve_interpolations && _interpolate_viscosity)
1275 19 : return _property_ve_ipol[_viscosity_idx]->sample(v, e);
1276 834 : else if (_construct_pT_from_ve)
1277 : {
1278 834 : Real p = _p_from_v_e_ipol->sample(v, e);
1279 834 : Real T = _T_from_v_e_ipol->sample(v, e);
1280 834 : return mu_from_p_T(p, T);
1281 : }
1282 0 : else if (_fp)
1283 0 : return _fp->mu_from_v_e(v, e);
1284 : else
1285 0 : mooseError(__PRETTY_FUNCTION__,
1286 : "\nNo tabulation or fluid property 'fp' object to compute value");
1287 : }
1288 :
1289 : void
1290 6 : TabulatedFluidProperties::mu_from_v_e(Real v, Real e, Real & mu, Real & dmu_dv, Real & dmu_de) const
1291 : {
1292 6 : if (_interpolate_viscosity & !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1293 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1294 6 : checkInputVariablesVE(v, e);
1295 :
1296 6 : if (_create_direct_ve_interpolations && _interpolate_viscosity)
1297 3 : _property_ve_ipol[_viscosity_idx]->sampleValueAndDerivatives(v, e, mu, dmu_dv, dmu_de);
1298 3 : else if (_construct_pT_from_ve)
1299 : {
1300 : Real p, dp_dv, dp_de;
1301 3 : _p_from_v_e_ipol->sampleValueAndDerivatives(v, e, p, dp_dv, dp_de);
1302 : Real T, dT_dv, dT_de;
1303 3 : _T_from_v_e_ipol->sampleValueAndDerivatives(v, e, T, dT_dv, dT_de);
1304 : Real dmu_dp, dmu_dT;
1305 3 : mu_from_p_T(p, T, mu, dmu_dp, dmu_dT);
1306 3 : dmu_dv = dmu_dp * dp_dv + dmu_dT * dT_dv;
1307 3 : dmu_de = dmu_dp * dp_de + dmu_dT * dT_de;
1308 : }
1309 0 : else if (_fp)
1310 0 : _fp->mu_from_v_e(v, e, mu, dmu_dv, dmu_de);
1311 : else
1312 0 : mooseError(__PRETTY_FUNCTION__,
1313 : "\nNo tabulation or fluid property 'fp' object to compute value");
1314 6 : }
1315 :
1316 : Real
1317 853 : TabulatedFluidProperties::k_from_v_e(Real v, Real e) const
1318 : {
1319 853 : if (_interpolate_k && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1320 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1321 853 : checkInputVariablesVE(v, e);
1322 :
1323 853 : if (_create_direct_ve_interpolations && _interpolate_k)
1324 19 : return _property_ve_ipol[_k_idx]->sample(v, e);
1325 834 : else if (_construct_pT_from_ve)
1326 : {
1327 834 : Real T = _T_from_v_e_ipol->sample(v, e);
1328 834 : Real p = _p_from_v_e_ipol->sample(v, e);
1329 834 : return k_from_p_T(p, T);
1330 : }
1331 0 : else if (_fp)
1332 0 : return _fp->k_from_v_e(v, e);
1333 : else
1334 0 : mooseError(__PRETTY_FUNCTION__,
1335 : "\nNo tabulation or fluid property 'fp' object to compute value");
1336 : }
1337 :
1338 : void
1339 6 : TabulatedFluidProperties::k_from_v_e(Real v, Real e, Real & k, Real & dk_dv, Real & dk_de) const
1340 : {
1341 6 : if (_interpolate_k && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1342 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1343 6 : checkInputVariablesVE(v, e);
1344 :
1345 6 : if (_create_direct_ve_interpolations && _interpolate_k)
1346 3 : _property_ve_ipol[_k_idx]->sampleValueAndDerivatives(v, e, k, dk_dv, dk_de);
1347 3 : else if (_construct_pT_from_ve)
1348 : {
1349 : Real p, dp_dv, dp_de;
1350 3 : _p_from_v_e_ipol->sampleValueAndDerivatives(v, e, p, dp_dv, dp_de);
1351 : Real T, dT_dv, dT_de;
1352 3 : _T_from_v_e_ipol->sampleValueAndDerivatives(v, e, T, dT_dv, dT_de);
1353 : Real dk_dp, dk_dT;
1354 3 : k_from_p_T(p, T, k, dk_dp, dk_dT);
1355 3 : dk_dv = dk_dp * dp_dv + dk_dT * dT_dv;
1356 3 : dk_de = dk_dp * dp_de + dk_dT * dT_de;
1357 : }
1358 0 : else if (_fp)
1359 0 : _fp->k_from_v_e(v, e, k, dk_dv, dk_de);
1360 : else
1361 0 : mooseError(__PRETTY_FUNCTION__,
1362 : "\nNo tabulation or fluid property 'fp' object to compute value");
1363 6 : }
1364 :
1365 : Real
1366 833 : TabulatedFluidProperties::s_from_v_e(Real v, Real e) const
1367 : {
1368 833 : if (_interpolate_entropy && !_construct_pT_from_ve && !_create_direct_ve_interpolations)
1369 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1370 833 : checkInputVariablesVE(v, e);
1371 :
1372 833 : if (_create_direct_ve_interpolations && _interpolate_entropy)
1373 9 : return _property_ve_ipol[_entropy_idx]->sample(v, e);
1374 824 : else if (_construct_pT_from_ve)
1375 : {
1376 824 : Real T = _T_from_v_e_ipol->sample(v, e);
1377 824 : Real p = _p_from_v_e_ipol->sample(v, e);
1378 824 : return s_from_p_T(p, T);
1379 : }
1380 0 : else if (_fp)
1381 0 : return _fp->s_from_v_e(v, e);
1382 : else
1383 0 : mooseError(__PRETTY_FUNCTION__,
1384 : "\nNo tabulation or fluid property 'fp' object to compute value");
1385 : }
1386 :
1387 : Real
1388 33 : TabulatedFluidProperties::g_from_v_e(Real v, Real e) const
1389 : {
1390 33 : if (!_construct_pT_from_ve &&
1391 9 : (!_create_direct_ve_interpolations || _entropy_idx == libMesh::invalid_uint ||
1392 9 : _enthalpy_idx == libMesh::invalid_uint || _T_idx == libMesh::invalid_uint))
1393 0 : missingVEInterpolationError(__PRETTY_FUNCTION__);
1394 33 : checkInputVariablesVE(v, e);
1395 :
1396 33 : Real h, T = 0, s;
1397 33 : if (_create_direct_ve_interpolations)
1398 : {
1399 9 : s = _property_ve_ipol[_entropy_idx]->sample(v, e);
1400 9 : h = _property_ve_ipol[_enthalpy_idx]->sample(v, e);
1401 9 : T = _property_ve_ipol[_T_idx]->sample(v, e);
1402 : }
1403 24 : else if (_fp || _create_direct_pT_interpolations)
1404 : {
1405 24 : Real p0 = _p_initial_guess;
1406 24 : Real T0 = _T_initial_guess;
1407 24 : Real p = 0;
1408 : bool conversion_succeeded;
1409 24 : p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded);
1410 24 : s = s_from_p_T(p, T);
1411 24 : h = h_from_p_T(p, T);
1412 24 : }
1413 : else
1414 0 : mooseError(__PRETTY_FUNCTION__,
1415 : "\nNo tabulation or fluid property 'fp' object to compute value");
1416 33 : return h - T * s;
1417 : }
1418 :
1419 : Real
1420 0 : TabulatedFluidProperties::T_from_h_s(Real h, Real s) const
1421 : {
1422 0 : Real p0 = _p_initial_guess;
1423 0 : Real T0 = _T_initial_guess;
1424 : Real p, T;
1425 : bool conversion_succeeded;
1426 0 : p_T_from_h_s(h, s, p0, T0, p, T, conversion_succeeded);
1427 0 : return T;
1428 : }
1429 :
1430 : Real
1431 4 : TabulatedFluidProperties::T_from_p_h(Real pressure, Real enthalpy) const
1432 : {
1433 4 : if (_fp)
1434 0 : return _fp->T_from_p_h(pressure, enthalpy);
1435 : else
1436 : {
1437 : auto lambda = [&](Real pressure, Real current_T, Real & new_h, Real & dh_dp, Real & dh_dT)
1438 9 : { h_from_p_T(pressure, current_T, new_h, dh_dp, dh_dT); };
1439 4 : Real T = FluidPropertiesUtils::NewtonSolve(
1440 4 : pressure, enthalpy, _T_initial_guess, _tolerance, lambda, name() + "::T_from_p_h")
1441 4 : .first;
1442 : // check for nans
1443 4 : if (std::isnan(T))
1444 0 : mooseError("Conversion from enthalpy (h = ",
1445 : enthalpy,
1446 : ") and pressure (p = ",
1447 : pressure,
1448 : ") to temperature failed to converge.");
1449 : return T;
1450 : }
1451 : }
1452 :
1453 : ADReal
1454 1 : TabulatedFluidProperties::T_from_p_h(const ADReal & pressure, const ADReal & enthalpy) const
1455 : {
1456 : using std::isnan;
1457 1 : if (_fp)
1458 0 : return _fp->T_from_p_h(pressure, enthalpy);
1459 : else
1460 : {
1461 : auto lambda =
1462 1 : [&](ADReal pressure, ADReal current_T, ADReal & new_h, ADReal & dh_dp, ADReal & dh_dT)
1463 : {
1464 1 : h_from_p_T(pressure.value(), current_T.value(), new_h.value(), dh_dp.value(), dh_dT.value());
1465 : // Reconstruct derivatives
1466 : new_h.derivatives() =
1467 1 : dh_dp.value() * pressure.derivatives() + dh_dT.value() * current_T.derivatives();
1468 1 : };
1469 : ADReal T =
1470 2 : FluidPropertiesUtils::NewtonSolve(
1471 2 : pressure, enthalpy, _T_initial_guess, _tolerance, lambda, name() + "::T_from_p_h")
1472 1 : .first;
1473 : // check for nans
1474 1 : if (isnan(T))
1475 0 : mooseError("Conversion from enthalpy (h = ",
1476 : enthalpy,
1477 : ") and pressure (p = ",
1478 : pressure,
1479 : ") to temperature failed to converge.");
1480 1 : return T;
1481 : }
1482 : }
1483 :
1484 : Real
1485 0 : TabulatedFluidProperties::s_from_h_p(Real enthalpy, Real pressure) const
1486 : {
1487 0 : Real T = T_from_p_h(pressure, enthalpy);
1488 0 : return s_from_p_T(pressure, T);
1489 : }
1490 :
1491 : void
1492 0 : TabulatedFluidProperties::s_from_h_p(
1493 : Real h, Real pressure, Real & s, Real & ds_dh, Real & ds_dp) const
1494 : {
1495 :
1496 0 : if (_fp)
1497 0 : _fp->s_from_h_p(h, pressure, s, ds_dh, ds_dp);
1498 : else
1499 0 : mooseError("fp", "s_from_h_p derivatives not implemented.");
1500 0 : }
1501 :
1502 : [[noreturn]] void
1503 0 : TabulatedFluidProperties::FluidPropertiesForwardError(const std::string & desired_routine) const
1504 : {
1505 0 : mooseError("TabulatedFluidProperties can only call the function '" + desired_routine +
1506 : "' when the 'fp' parameter is provided. It is currently not implemented using "
1507 : "tabulations, and this property is simply forwarded to the FluidProperties specified "
1508 : "in the 'fp' parameter");
1509 : }
1510 :
1511 : void
1512 24 : TabulatedFluidProperties::writeTabulatedData(std::string file_name)
1513 : {
1514 48 : file_name = file_name.empty() ? "fluid_properties_" + name() + "_out.csv" : file_name;
1515 24 : if (processor_id() == 0)
1516 : {
1517 : {
1518 16 : MooseUtils::checkFileWriteable(file_name);
1519 :
1520 16 : std::ofstream file_out(file_name.c_str());
1521 :
1522 : // Write out date and fluid type
1523 16 : time_t now = std::time(&now);
1524 16 : if (_fp)
1525 16 : file_out << "# " << _fp->fluidName()
1526 16 : << " properties created by TabulatedFluidProperties on " << ctime(&now) << "\n";
1527 : else
1528 8 : file_out << "# tabulated properties created by TabulatedFluidProperties on " << ctime(&now)
1529 16 : << "\n";
1530 :
1531 : // Write out column names
1532 16 : file_out << "pressure, temperature";
1533 160 : for (std::size_t i = 0; i < _interpolated_properties.size(); ++i)
1534 144 : file_out << ", " << _interpolated_properties[i];
1535 16 : file_out << "\n";
1536 :
1537 : // Write out the fluid property data
1538 1616 : for (unsigned int p = 0; p < _num_p; ++p)
1539 161600 : for (unsigned int t = 0; t < _num_T; ++t)
1540 : {
1541 160000 : file_out << _pressure[p] << ", " << _temperature[t];
1542 1600000 : for (std::size_t i = 0; i < _properties.size(); ++i)
1543 1440000 : file_out << ", " << _properties[i][p * _num_T + t];
1544 160000 : file_out << "\n";
1545 : }
1546 :
1547 : file_out << std::flush;
1548 16 : file_out.close();
1549 16 : }
1550 :
1551 : // Write out the (v,e) to (p,T) conversions
1552 16 : if (_construct_pT_from_ve)
1553 : {
1554 8 : const auto file_name_ve = (_file_name_ve_out == "")
1555 16 : ? std::regex_replace(file_name, std::regex("\\.csv"), "_ve.csv")
1556 8 : : _file_name_ve_out;
1557 8 : MooseUtils::checkFileWriteable(file_name_ve);
1558 8 : std::ofstream file_out(file_name_ve.c_str());
1559 :
1560 : // Write out column names
1561 8 : file_out << "specific_volume, internal_energy, pressure, temperature";
1562 80 : for (const auto i : index_range(_properties))
1563 72 : if (_interpolated_properties[i] != "internal_energy")
1564 64 : file_out << ", " << _interpolated_properties[i];
1565 8 : file_out << "\n";
1566 :
1567 : // Write out the fluid property data
1568 88 : for (const auto v : make_range(_num_v))
1569 880 : for (const auto e : make_range(_num_e))
1570 : {
1571 800 : const auto v_val = _specific_volume[v];
1572 800 : const auto e_val = _internal_energy[e];
1573 800 : const auto pressure = _p_from_v_e_ipol->sample(v_val, e_val);
1574 800 : const auto temperature = _T_from_v_e_ipol->sample(v_val, e_val);
1575 3200 : file_out << v_val << ", " << e_val << ", " << pressure << ", " << temperature << ", ";
1576 8000 : for (const auto i : index_range(_properties))
1577 : {
1578 : bool add_comma = true;
1579 7200 : if (i == _density_idx)
1580 800 : file_out << 1 / v_val;
1581 6400 : else if (i == _enthalpy_idx)
1582 800 : file_out << h_from_p_T(pressure, temperature);
1583 : // Note that we could use (p,T) routine to generate this instead of (v,e)
1584 : // Or could use the _properties_ve array
1585 5600 : else if (i == _viscosity_idx)
1586 800 : file_out << mu_from_v_e(v_val, e_val);
1587 4800 : else if (i == _k_idx)
1588 800 : file_out << k_from_v_e(v_val, e_val);
1589 4000 : else if (i == _c_idx)
1590 800 : file_out << c_from_v_e(v_val, e_val);
1591 3200 : else if (i == _cv_idx)
1592 800 : file_out << cv_from_v_e(v_val, e_val);
1593 2400 : else if (i == _cp_idx)
1594 800 : file_out << cp_from_v_e(v_val, e_val);
1595 1600 : else if (i == _entropy_idx)
1596 800 : file_out << s_from_v_e(v_val, e_val);
1597 : else
1598 : add_comma = false;
1599 7200 : if (i != _properties.size() - 1 && add_comma)
1600 5600 : file_out << ", ";
1601 : }
1602 :
1603 800 : file_out << "\n";
1604 : }
1605 :
1606 : file_out << std::flush;
1607 8 : file_out.close();
1608 8 : }
1609 : }
1610 24 : }
1611 :
1612 : void
1613 126 : TabulatedFluidProperties::generateTabulatedData()
1614 : {
1615 : mooseAssert(_fp, "We should not try to generate (p,T) tabulated data without a _fp user object");
1616 126 : _pressure.resize(_num_p);
1617 126 : _temperature.resize(_num_T);
1618 :
1619 : // Generate data for all properties entered in input file
1620 126 : _properties.resize(_interpolated_properties_enum.size());
1621 126 : _interpolated_properties.resize(_interpolated_properties_enum.size());
1622 :
1623 1144 : for (std::size_t i = 0; i < _interpolated_properties_enum.size(); ++i)
1624 1018 : _interpolated_properties[i] = _interpolated_properties_enum[i];
1625 :
1626 1144 : for (const auto i : index_range(_properties))
1627 1018 : _properties[i].resize(_num_p * _num_T);
1628 :
1629 : // Temperature is divided equally into _num_T segments
1630 126 : Real delta_T = (_temperature_max - _temperature_min) / static_cast<Real>(_num_T - 1);
1631 :
1632 12632 : for (unsigned int j = 0; j < _num_T; ++j)
1633 12506 : _temperature[j] = _temperature_min + j * delta_T;
1634 :
1635 : // Divide the pressure into _num_p equal segments
1636 126 : Real delta_p = (_pressure_max - _pressure_min) / static_cast<Real>(_num_p - 1);
1637 :
1638 12632 : for (unsigned int i = 0; i < _num_p; ++i)
1639 12506 : _pressure[i] = _pressure_min + i * delta_p;
1640 :
1641 : // Generate the tabulated data at the pressure and temperature points
1642 1144 : for (const auto i : index_range(_properties))
1643 : {
1644 1018 : if (_interpolated_properties[i] == "density")
1645 12632 : for (unsigned int p = 0; p < _num_p; ++p)
1646 1262542 : for (unsigned int t = 0; t < _num_T; ++t)
1647 1250036 : _properties[i][p * _num_T + t] = _fp->rho_from_p_T(_pressure[p], _temperature[t]);
1648 :
1649 1018 : if (_interpolated_properties[i] == "enthalpy")
1650 12632 : for (unsigned int p = 0; p < _num_p; ++p)
1651 1262542 : for (unsigned int t = 0; t < _num_T; ++t)
1652 1250036 : _properties[i][p * _num_T + t] = _fp->h_from_p_T(_pressure[p], _temperature[t]);
1653 :
1654 1018 : if (_interpolated_properties[i] == "internal_energy")
1655 12632 : for (unsigned int p = 0; p < _num_p; ++p)
1656 1262542 : for (unsigned int t = 0; t < _num_T; ++t)
1657 1250036 : _properties[i][p * _num_T + t] = _fp->e_from_p_T(_pressure[p], _temperature[t]);
1658 :
1659 1018 : if (_interpolated_properties[i] == "viscosity")
1660 12632 : for (unsigned int p = 0; p < _num_p; ++p)
1661 1262542 : for (unsigned int t = 0; t < _num_T; ++t)
1662 1250036 : _properties[i][p * _num_T + t] = _fp->mu_from_p_T(_pressure[p], _temperature[t]);
1663 :
1664 1018 : if (_interpolated_properties[i] == "k")
1665 10309 : for (unsigned int p = 0; p < _num_p; ++p)
1666 1030242 : for (unsigned int t = 0; t < _num_T; ++t)
1667 1020036 : _properties[i][p * _num_T + t] = _fp->k_from_p_T(_pressure[p], _temperature[t]);
1668 :
1669 1018 : if (_interpolated_properties[i] == "c")
1670 10302 : for (unsigned int p = 0; p < _num_p; ++p)
1671 1030200 : for (unsigned int t = 0; t < _num_T; ++t)
1672 1020000 : _properties[i][p * _num_T + t] = _fp->c_from_p_T(_pressure[p], _temperature[t]);
1673 :
1674 1018 : if (_interpolated_properties[i] == "cv")
1675 10309 : for (unsigned int p = 0; p < _num_p; ++p)
1676 1030242 : for (unsigned int t = 0; t < _num_T; ++t)
1677 1020036 : _properties[i][p * _num_T + t] = _fp->cv_from_p_T(_pressure[p], _temperature[t]);
1678 :
1679 1018 : if (_interpolated_properties[i] == "cp")
1680 10309 : for (unsigned int p = 0; p < _num_p; ++p)
1681 1030242 : for (unsigned int t = 0; t < _num_T; ++t)
1682 1020036 : _properties[i][p * _num_T + t] = _fp->cp_from_p_T(_pressure[p], _temperature[t]);
1683 :
1684 1018 : if (_interpolated_properties[i] == "entropy")
1685 10309 : for (unsigned int p = 0; p < _num_p; ++p)
1686 1030242 : for (unsigned int t = 0; t < _num_T; ++t)
1687 1020036 : _properties[i][p * _num_T + t] = _fp->s_from_p_T(_pressure[p], _temperature[t]);
1688 : }
1689 126 : }
1690 :
1691 : void
1692 2 : TabulatedFluidProperties::generateVETabulatedData()
1693 : {
1694 : mooseAssert(_fp, "We should not try to generate (v,e) tabulated data without a _fp user object");
1695 2 : _specific_volume.resize(_num_v);
1696 2 : _internal_energy.resize(_num_e);
1697 :
1698 : // Generate data for all properties entered in input file
1699 2 : _properties_ve.resize(_interpolated_properties_enum.size());
1700 2 : _interpolated_properties.resize(_interpolated_properties_enum.size());
1701 :
1702 : // This is filled from the user input, so it does not matter than this operation is performed
1703 : // for both the (p,T) and (v,e) tabulated data generation
1704 22 : for (std::size_t i = 0; i < _interpolated_properties_enum.size(); ++i)
1705 20 : _interpolated_properties[i] = _interpolated_properties_enum[i];
1706 :
1707 22 : for (const auto i : index_range(_properties_ve))
1708 20 : _properties_ve[i].resize(_num_v * _num_e);
1709 :
1710 : // Grids in (v,e) are not read, so we either use user input or rely on (p,T) data
1711 2 : createVEGridVectors();
1712 :
1713 : // Generate the tabulated data at the pressure and temperature points
1714 22 : for (const auto i : index_range(_properties_ve))
1715 : {
1716 20 : if (_interpolated_properties[i] == "density")
1717 202 : for (unsigned int v = 0; v < _num_v; ++v)
1718 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1719 20000 : _properties_ve[i][v * _num_e + e] = 1. / _specific_volume[v];
1720 :
1721 20 : if (_interpolated_properties[i] == "enthalpy")
1722 202 : for (unsigned int v = 0; v < _num_v; ++v)
1723 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1724 20000 : _properties_ve[i][v * _num_e + e] =
1725 20000 : _fp->h_from_v_e(_specific_volume[v], _internal_energy[e]);
1726 :
1727 20 : if (_interpolated_properties[i] == "internal_energy")
1728 0 : for (unsigned int v = 0; v < _num_v; ++v)
1729 0 : for (unsigned int e = 0; e < _num_e; ++e)
1730 0 : _properties_ve[i][v * _num_e + e] = _internal_energy[e];
1731 :
1732 20 : if (_interpolated_properties[i] == "viscosity")
1733 202 : for (unsigned int v = 0; v < _num_v; ++v)
1734 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1735 20000 : _properties_ve[i][v * _num_e + e] =
1736 20000 : _fp->mu_from_v_e(_specific_volume[v], _internal_energy[e]);
1737 :
1738 20 : if (_interpolated_properties[i] == "k")
1739 202 : for (unsigned int v = 0; v < _num_v; ++v)
1740 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1741 20000 : _properties_ve[i][v * _num_e + e] =
1742 20000 : _fp->k_from_v_e(_specific_volume[v], _internal_energy[e]);
1743 :
1744 20 : if (_interpolated_properties[i] == "c")
1745 202 : for (unsigned int v = 0; v < _num_v; ++v)
1746 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1747 20000 : _properties_ve[i][v * _num_e + e] =
1748 20000 : _fp->c_from_v_e(_specific_volume[v], _internal_energy[e]);
1749 :
1750 20 : if (_interpolated_properties[i] == "cv")
1751 202 : for (unsigned int v = 0; v < _num_v; ++v)
1752 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1753 20000 : _properties_ve[i][v * _num_e + e] =
1754 20000 : _fp->cv_from_v_e(_specific_volume[v], _internal_energy[e]);
1755 :
1756 20 : if (_interpolated_properties[i] == "cp")
1757 202 : for (unsigned int v = 0; v < _num_v; ++v)
1758 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1759 20000 : _properties_ve[i][v * _num_e + e] =
1760 20000 : _fp->cp_from_v_e(_specific_volume[v], _internal_energy[e]);
1761 :
1762 20 : if (_interpolated_properties[i] == "entropy")
1763 202 : for (unsigned int v = 0; v < _num_v; ++v)
1764 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1765 20000 : _properties_ve[i][v * _num_e + e] =
1766 20000 : _fp->s_from_v_e(_specific_volume[v], _internal_energy[e]);
1767 :
1768 20 : if (_interpolated_properties[i] == "pressure")
1769 202 : for (unsigned int v = 0; v < _num_v; ++v)
1770 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1771 20000 : _properties_ve[i][v * _num_e + e] =
1772 20000 : _fp->p_from_v_e(_specific_volume[v], _internal_energy[e]);
1773 :
1774 20 : if (_interpolated_properties[i] == "temperature")
1775 202 : for (unsigned int v = 0; v < _num_v; ++v)
1776 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1777 20000 : _properties_ve[i][v * _num_e + e] =
1778 20000 : _fp->T_from_v_e(_specific_volume[v], _internal_energy[e]);
1779 : }
1780 2 : }
1781 :
1782 : template <typename T>
1783 : void
1784 37000049 : TabulatedFluidProperties::checkInputVariables(T & pressure, T & temperature) const
1785 : {
1786 : using std::max, std::min;
1787 :
1788 37000049 : if (_OOBBehavior == Ignore)
1789 : return;
1790 36999401 : else if (MooseUtils::absoluteFuzzyGreaterThan(_pressure_min, pressure, libMesh::TOLERANCE) ||
1791 : MooseUtils::absoluteFuzzyGreaterThan(pressure, _pressure_max, libMesh::TOLERANCE))
1792 : {
1793 28249314 : if (_OOBBehavior == Throw)
1794 12 : throw MooseException("Pressure " + Moose::stringify(pressure) +
1795 : " is outside the range of tabulated pressure (" +
1796 2 : Moose::stringify(_pressure_min) + ", " +
1797 2 : Moose::stringify(_pressure_max) + ").");
1798 :
1799 : else
1800 : {
1801 36915376 : pressure = max(_pressure_min, min(pressure, _pressure_max));
1802 28249312 : if (_OOBBehavior == DeclareInvalid)
1803 146 : flagInvalidSolution("Pressure out of bounds");
1804 28249168 : else if (_OOBBehavior == WarnInvalid)
1805 545 : flagSolutionWarning("Pressure out of bounds");
1806 : }
1807 : }
1808 :
1809 36999399 : if (MooseUtils::absoluteFuzzyGreaterThan(_temperature_min, temperature, libMesh::TOLERANCE) ||
1810 : MooseUtils::absoluteFuzzyGreaterThan(temperature, _temperature_max, libMesh::TOLERANCE))
1811 : {
1812 11054558 : if (_OOBBehavior == Throw)
1813 12 : throw MooseException("Temperature " + Moose::stringify(temperature) +
1814 : " is outside the range of tabulated temperature (" +
1815 2 : Moose::stringify(_temperature_min) + ", " +
1816 2 : Moose::stringify(_temperature_max) + ").");
1817 : else
1818 : {
1819 12574348 : temperature = max(T(_temperature_min), min(temperature, T(_temperature_max)));
1820 11054556 : if (_OOBBehavior == DeclareInvalid)
1821 146 : flagInvalidSolution("Temperature out of bounds");
1822 11054412 : else if (_OOBBehavior == WarnInvalid)
1823 545 : flagSolutionWarning("Temperature out of bounds");
1824 : }
1825 : }
1826 : }
1827 :
1828 : template <typename T>
1829 : void
1830 5285 : TabulatedFluidProperties::checkInputVariablesVE(T & v, T & e) const
1831 : {
1832 : using std::max, std::min;
1833 :
1834 5285 : if (_OOBBehavior == Ignore)
1835 : return;
1836 5285 : else if (e < _e_min || e > _e_max)
1837 : {
1838 67 : if (_OOBBehavior == Throw)
1839 0 : throw MooseException("Specific internal energy " + Moose::stringify(e) +
1840 : " is outside the range of tabulated specific internal energies (" +
1841 0 : Moose::stringify(_e_min) + ", " + Moose::stringify(_e_max) + ").");
1842 : else
1843 : {
1844 80 : e = max(_e_min, min(e, _e_max));
1845 67 : if (_OOBBehavior == DeclareInvalid)
1846 0 : flagInvalidSolution("Specific internal energy out of bounds");
1847 67 : else if (_OOBBehavior == WarnInvalid)
1848 0 : flagSolutionWarning("Specific internal energy out of bounds");
1849 : }
1850 : }
1851 :
1852 5285 : if (v < _v_min || v > _v_max)
1853 : {
1854 54 : if (_OOBBehavior == Throw)
1855 0 : throw MooseException("Specific volume " + Moose::stringify(v) +
1856 : " is outside the range of tabulated specific volumes (" +
1857 0 : Moose::stringify(_v_min) + ", " + Moose::stringify(_v_max) + ").");
1858 : else
1859 : {
1860 54 : v = max(T(_v_min), min(v, T(_v_max)));
1861 54 : if (_OOBBehavior == DeclareInvalid)
1862 0 : flagInvalidSolution("Specific volume out of bounds");
1863 54 : else if (_OOBBehavior == WarnInvalid)
1864 0 : flagSolutionWarning("Specific volume out of bounds");
1865 : }
1866 : }
1867 : }
1868 :
1869 : void
1870 326 : TabulatedFluidProperties::checkInitialGuess() const
1871 : {
1872 326 : if (_fp && (_construct_pT_from_ve || _construct_pT_from_vh))
1873 : {
1874 16 : if (_p_initial_guess < _pressure_min || _p_initial_guess > _pressure_max)
1875 28 : mooseWarning("Pressure initial guess for (p,T), (v,e) conversions " +
1876 26 : Moose::stringify(_p_initial_guess) +
1877 : " is outside the range of tabulated "
1878 14 : "pressure (" +
1879 40 : Moose::stringify(_pressure_min) + ", " + Moose::stringify(_pressure_max) + ").");
1880 :
1881 14 : if (_T_initial_guess < _temperature_min || _T_initial_guess > _temperature_max)
1882 28 : mooseWarning("Temperature initial guess for (p,T), (v,e) conversions " +
1883 26 : Moose::stringify(_T_initial_guess) +
1884 : " is outside the range of tabulated "
1885 14 : "temperature (" +
1886 40 : Moose::stringify(_temperature_min) + ", " + Moose::stringify(_temperature_max) +
1887 : ").");
1888 : }
1889 322 : }
1890 :
1891 : void
1892 65 : TabulatedFluidProperties::readFileTabulationData(const bool use_pT)
1893 : {
1894 : std::string file_name;
1895 65 : if (use_pT)
1896 : {
1897 43 : _console << name() + ": Reading tabulated properties from " << _file_name_in << std::endl;
1898 43 : _csv_reader.read();
1899 43 : file_name = _file_name_in;
1900 : }
1901 : else
1902 : {
1903 27 : _console << name() + ": Reading tabulated properties from " << _file_name_ve_in << std::endl;
1904 22 : _csv_reader.setFileName(_file_name_ve_in);
1905 22 : _csv_reader.read();
1906 : file_name = _file_name_ve_in;
1907 : }
1908 :
1909 65 : const std::vector<std::string> & column_names = _csv_reader.getNames();
1910 :
1911 : // These columns form the grid and must be present in the file
1912 : std::vector<std::string> required_columns;
1913 65 : if (use_pT)
1914 129 : required_columns = {"pressure", "temperature"};
1915 : else
1916 66 : required_columns = {"specific_volume", "internal_energy"};
1917 :
1918 : // Check that all required columns are present
1919 194 : for (std::size_t i = 0; i < required_columns.size(); ++i)
1920 : {
1921 130 : if (std::find(column_names.begin(), column_names.end(), required_columns[i]) ==
1922 : column_names.end())
1923 1 : mooseError("No ",
1924 : required_columns[i],
1925 : " data read in ",
1926 : file_name,
1927 : ". A column named ",
1928 : required_columns[i],
1929 : " must be present");
1930 : }
1931 :
1932 : // These columns can be present in the file
1933 : std::vector<std::string> property_columns = {
1934 64 : "density", "enthalpy", "viscosity", "k", "c", "cv", "cp", "entropy"};
1935 64 : if (use_pT)
1936 84 : property_columns.push_back("internal_energy");
1937 : else
1938 : {
1939 22 : property_columns.push_back("pressure");
1940 44 : property_columns.push_back("temperature");
1941 : }
1942 :
1943 : // Check that any property names read from the file are present in the list of possible
1944 : // properties, and if they are, add them to the list of read properties
1945 762 : for (std::size_t i = 0; i < column_names.size(); ++i)
1946 : {
1947 : // Only check properties not in _required_columns
1948 699 : if (std::find(required_columns.begin(), required_columns.end(), column_names[i]) ==
1949 : required_columns.end())
1950 : {
1951 571 : if (std::find(property_columns.begin(), property_columns.end(), column_names[i]) ==
1952 : property_columns.end())
1953 1 : mooseWarning(column_names[i],
1954 : " read in ",
1955 : file_name,
1956 : " tabulation file is not one of the properties that TabulatedFluidProperties "
1957 : "understands. It will be ignored.");
1958 : // We could be reading a (v,e) tabulation after having read a (p,T) tabulation, do not
1959 : // insert twice
1960 570 : else if (std::find(_interpolated_properties.begin(),
1961 : _interpolated_properties.end(),
1962 : column_names[i]) == _interpolated_properties.end())
1963 570 : _interpolated_properties.push_back(column_names[i]);
1964 : }
1965 : }
1966 :
1967 : std::map<std::string, unsigned int> data_index;
1968 758 : for (std::size_t i = 0; i < column_names.size(); ++i)
1969 : {
1970 695 : auto it = std::find(column_names.begin(), column_names.end(), column_names[i]);
1971 695 : data_index[column_names[i]] = std::distance(column_names.begin(), it);
1972 : }
1973 :
1974 63 : const std::vector<std::vector<Real>> & column_data = _csv_reader.getData();
1975 :
1976 : // Extract the pressure and temperature data vectors
1977 63 : if (use_pT)
1978 : {
1979 41 : _pressure = column_data[data_index.find("pressure")->second];
1980 82 : _temperature = column_data[data_index.find("temperature")->second];
1981 : }
1982 : else
1983 : {
1984 22 : _specific_volume = column_data[data_index.find("specific_volume")->second];
1985 44 : _internal_energy = column_data[data_index.find("internal_energy")->second];
1986 : }
1987 :
1988 63 : if (use_pT)
1989 82 : checkFileTabulationGrids(_pressure, _temperature, file_name, "pressure", "temperature");
1990 : else
1991 44 : checkFileTabulationGrids(_specific_volume,
1992 22 : _internal_energy,
1993 : file_name,
1994 : "specific volume",
1995 : "specific internal energy");
1996 :
1997 60 : if (use_pT)
1998 : {
1999 38 : _num_p = _pressure.size();
2000 38 : _num_T = _temperature.size();
2001 :
2002 : // Minimum and maximum pressure and temperature. Note that _pressure and
2003 : // _temperature are sorted
2004 38 : _pressure_min = _pressure.front();
2005 38 : _pressure_max = _pressure.back();
2006 38 : _temperature_min = _temperature.front();
2007 38 : _temperature_max = _temperature.back();
2008 :
2009 : // Extract the fluid property data from the file
2010 378 : for (std::size_t i = 0; i < _interpolated_properties.size(); ++i)
2011 340 : _properties.push_back(column_data[data_index.find(_interpolated_properties[i])->second]);
2012 : }
2013 : else
2014 : {
2015 22 : _num_v = _specific_volume.size();
2016 22 : _num_e = _internal_energy.size();
2017 :
2018 : // Minimum and maximum specific internal energy and specific volume
2019 22 : _v_min = _specific_volume.front();
2020 22 : _v_max = _specific_volume.back();
2021 22 : _e_min = _internal_energy.front();
2022 22 : _e_max = _internal_energy.back();
2023 :
2024 : // We cannot overwrite the tabulated data grid with a grid generated from user-input for the
2025 : // purpose of creating (p,T) to (v,e) interpolations
2026 22 : if (_construct_pT_from_ve)
2027 3 : paramError("construct_pT_from_ve",
2028 : "Reading a (v,e) tabulation and generating (p,T) to (v,e) interpolation tables is "
2029 : "not supported at this time.");
2030 :
2031 : // Make sure we use the tabulation bounds
2032 22 : _e_bounds_specified = true;
2033 22 : _v_bounds_specified = true;
2034 :
2035 : // Extract the fluid property data from the file
2036 22 : _properties_ve.reserve(_interpolated_properties.size());
2037 242 : for (std::size_t i = 0; i < _interpolated_properties.size(); ++i)
2038 220 : _properties_ve.push_back(column_data[data_index.find(_interpolated_properties[i])->second]);
2039 : }
2040 259 : }
2041 :
2042 : void
2043 63 : TabulatedFluidProperties::checkFileTabulationGrids(std::vector<Real> & v1,
2044 : std::vector<Real> & v2,
2045 : const std::string & file_name,
2046 : const std::string & v1_name,
2047 : const std::string & v2_name)
2048 : {
2049 : // NOTE: We kept the comments in terms of pressure & temperature for clarity
2050 : // Pressure (v1) and temperature (v2) data contains duplicates due to the csv format.
2051 : // First, check that pressure (v1) is monotonically increasing
2052 63 : if (!std::is_sorted(v1.begin(), v1.end()))
2053 0 : mooseError("The column data for ", v1_name, " is not monotonically increasing in ", file_name);
2054 :
2055 : // The first pressure (v1) value is repeated for each temperature (v2) value. Counting the
2056 : // number of repeats provides the number of temperature (v2) values
2057 63 : auto num_v2 = std::count(v1.begin(), v1.end(), v1.front());
2058 :
2059 : // Now remove the duplicates in the pressure (v1) vector
2060 63 : auto last_unique = std::unique(v1.begin(), v1.end());
2061 : v1.erase(last_unique, v1.end());
2062 :
2063 : // Check that the number of rows in the csv file is equal to _num_v1 * _num_v2
2064 : // v2 is currently the same size as the column_data (will get trimmed at the end)
2065 63 : if (v2.size() != v1.size() * libMesh::cast_int<unsigned int>(num_v2))
2066 1 : mooseError("The number of rows in ",
2067 : file_name,
2068 : " is not equal to the number of unique ",
2069 : v1_name,
2070 : " values ",
2071 1 : v1.size(),
2072 : " multiplied by the number of unique ",
2073 : v2_name,
2074 : " values ",
2075 : num_v2);
2076 :
2077 : // Need to make sure that the temperature (v2) values are provided in ascending order
2078 62 : std::vector<Real> base_v2(v2.begin(), v2.begin() + num_v2);
2079 62 : if (!std::is_sorted(base_v2.begin(), base_v2.end()))
2080 1 : mooseError("The column data for ", v2_name, " is not monotonically increasing in ", file_name);
2081 :
2082 : // Need to make sure that the temperature (v2) are repeated for each pressure (v1) grid point
2083 61 : auto it_v2 = v2.begin() + num_v2;
2084 3832 : for (const auto i : make_range(v1.size() - 1))
2085 : {
2086 3773 : std::vector<Real> repeated_v2(it_v2, it_v2 + num_v2);
2087 : if (repeated_v2 != base_v2)
2088 1 : mooseError(v2_name,
2089 : " values for ",
2090 : v1_name,
2091 : " ",
2092 1 : v1[i + 1],
2093 : " are not identical to values for ",
2094 : v1[0]);
2095 :
2096 : std::advance(it_v2, num_v2);
2097 3772 : }
2098 :
2099 : // At this point, all temperature (v2) data has been provided in ascending order
2100 : // identically for each pressure (v1) value, so we can just keep the first range
2101 : v2.erase(v2.begin() + num_v2, v2.end());
2102 62 : }
2103 :
2104 : void
2105 188 : TabulatedFluidProperties::computePropertyIndicesInInterpolationVectors()
2106 : {
2107 : // At this point, all properties read or generated are able to be used by
2108 : // TabulatedFluidProperties. Now set flags and indexes for each property in
2109 : //_interpolated_properties to use in property calculations
2110 1786 : for (std::size_t i = 0; i < _interpolated_properties.size(); ++i)
2111 : {
2112 1598 : if (_interpolated_properties[i] == "density")
2113 : {
2114 188 : _interpolate_density = true;
2115 188 : _density_idx = i;
2116 : }
2117 1410 : else if (_interpolated_properties[i] == "enthalpy")
2118 : {
2119 188 : _interpolate_enthalpy = true;
2120 188 : _enthalpy_idx = i;
2121 : }
2122 1222 : else if (_interpolated_properties[i] == "internal_energy")
2123 : {
2124 164 : _interpolate_internal_energy = true;
2125 164 : _internal_energy_idx = i;
2126 : }
2127 1058 : else if (_interpolated_properties[i] == "viscosity")
2128 : {
2129 188 : _interpolate_viscosity = true;
2130 188 : _viscosity_idx = i;
2131 : }
2132 870 : else if (_interpolated_properties[i] == "k")
2133 : {
2134 165 : _interpolate_k = true;
2135 165 : _k_idx = i;
2136 : }
2137 705 : else if (_interpolated_properties[i] == "c")
2138 : {
2139 162 : _interpolate_c = true;
2140 162 : _c_idx = i;
2141 : }
2142 543 : else if (_interpolated_properties[i] == "cp")
2143 : {
2144 165 : _interpolate_cp = true;
2145 165 : _cp_idx = i;
2146 : }
2147 378 : else if (_interpolated_properties[i] == "cv")
2148 : {
2149 165 : _interpolate_cv = true;
2150 165 : _cv_idx = i;
2151 : }
2152 213 : else if (_interpolated_properties[i] == "entropy")
2153 : {
2154 165 : _interpolate_entropy = true;
2155 165 : _entropy_idx = i;
2156 : }
2157 48 : else if (_interpolated_properties[i] == "pressure")
2158 : {
2159 24 : _interpolate_pressure = true;
2160 24 : _p_idx = i;
2161 : }
2162 24 : else if (_interpolated_properties[i] == "temperature")
2163 : {
2164 24 : _interpolate_temperature = true;
2165 24 : _T_idx = i;
2166 : }
2167 : else
2168 0 : mooseError("Specified property '" + _interpolated_properties[i] +
2169 : "' is present in the tabulation but is not currently leveraged by the code in the "
2170 : "TabulatedFluidProperties. If it is spelled correctly, then please contact a "
2171 : "MOOSE or fluid properties module developer.");
2172 : }
2173 188 : }
2174 :
2175 : void
2176 102 : TabulatedFluidProperties::createVGridVector()
2177 : {
2178 : mooseAssert(_file_name_ve_in.empty(), "We should be reading the specific volume grid from file");
2179 102 : if (!_v_bounds_specified)
2180 : {
2181 52 : if (_fp)
2182 : {
2183 : // extreme values of specific volume for the grid bounds
2184 14 : Real v1 = v_from_p_T(_pressure_min, _temperature_min);
2185 14 : Real v2 = v_from_p_T(_pressure_max, _temperature_min);
2186 14 : Real v3 = v_from_p_T(_pressure_min, _temperature_max);
2187 14 : Real v4 = v_from_p_T(_pressure_max, _temperature_max);
2188 14 : _v_min = std::min({v1, v2, v3, v4});
2189 14 : _v_max = std::max({v1, v2, v3, v4});
2190 : }
2191 : // if csv exists, get max and min values from csv file
2192 : else
2193 : {
2194 : Real rho_max =
2195 76 : *std::max_element(_properties[_density_idx].begin(), _properties[_density_idx].end());
2196 : Real rho_min =
2197 38 : *std::min_element(_properties[_density_idx].begin(), _properties[_density_idx].end());
2198 38 : _v_max = 1 / rho_min;
2199 38 : _v_min = 1 / rho_max;
2200 : }
2201 : // Prevent changing the bounds of the grid
2202 52 : _v_bounds_specified = true;
2203 : }
2204 :
2205 : // Create v grid for interpolation
2206 102 : _specific_volume.resize(_num_v);
2207 102 : if (_log_space_v)
2208 : {
2209 : // incrementing the exponent linearly will yield a log-spaced grid after taking the value to
2210 : // the power of 10
2211 24 : Real dv = (std::log10(_v_max) - std::log10(_v_min)) / ((Real)_num_v - 1);
2212 24 : Real log_v_min = std::log10(_v_min);
2213 2424 : for (unsigned int j = 0; j < _num_v; ++j)
2214 2400 : _specific_volume[j] = std::pow(10, log_v_min + j * dv);
2215 : }
2216 : else
2217 : {
2218 78 : Real dv = (_v_max - _v_min) / ((Real)_num_v - 1);
2219 5718 : for (unsigned int j = 0; j < _num_v; ++j)
2220 5640 : _specific_volume[j] = _v_min + j * dv;
2221 : }
2222 102 : }
2223 :
2224 : void
2225 52 : TabulatedFluidProperties::createVEGridVectors()
2226 : {
2227 52 : createVGridVector();
2228 52 : if (!_e_bounds_specified)
2229 : {
2230 52 : if (_fp)
2231 : {
2232 : // extreme values of internal energy for the grid bounds
2233 14 : Real e1 = e_from_p_T(_pressure_min, _temperature_min);
2234 14 : Real e2 = e_from_p_T(_pressure_max, _temperature_min);
2235 14 : Real e3 = e_from_p_T(_pressure_min, _temperature_max);
2236 14 : Real e4 = e_from_p_T(_pressure_max, _temperature_max);
2237 14 : _e_min = std::min({e1, e2, e3, e4});
2238 14 : _e_max = std::max({e1, e2, e3, e4});
2239 : }
2240 : // if csv exists, get max and min values from csv file
2241 : else
2242 : {
2243 38 : _e_min = *std::min_element(_properties[_internal_energy_idx].begin(),
2244 38 : _properties[_internal_energy_idx].end());
2245 38 : _e_max = *std::max_element(_properties[_internal_energy_idx].begin(),
2246 : _properties[_internal_energy_idx].end());
2247 : }
2248 : }
2249 :
2250 : // Create e grid for interpolation
2251 52 : _internal_energy.resize(_num_e);
2252 52 : if (_log_space_e)
2253 : {
2254 : // incrementing the exponent linearly will yield a log-spaced grid after taking the value to
2255 : // the power of 10
2256 12 : if (_e_min < 0)
2257 0 : mooseError("Logarithmic grid in specific energy can only be used with a positive specific "
2258 0 : "energy. Current minimum: " +
2259 : std::to_string(_e_min));
2260 12 : Real de = (std::log10(_e_max) - std::log10(_e_min)) / ((Real)_num_e - 1);
2261 12 : Real log_e_min = std::log10(_e_min);
2262 1212 : for (const auto j : make_range(_num_e))
2263 1200 : _internal_energy[j] = std::pow(10, log_e_min + j * de);
2264 : }
2265 : else
2266 : {
2267 40 : Real de = (_e_max - _e_min) / ((Real)_num_e - 1);
2268 2960 : for (const auto j : make_range(_num_e))
2269 2920 : _internal_energy[j] = _e_min + j * de;
2270 : }
2271 52 : }
2272 :
2273 : void
2274 50 : TabulatedFluidProperties::createVHGridVectors()
2275 : {
2276 50 : if (_file_name_ve_in.empty())
2277 50 : createVGridVector();
2278 50 : if (_fp)
2279 : {
2280 : // extreme values of enthalpy for the grid bounds
2281 12 : Real h1 = h_from_p_T(_pressure_min, _temperature_min);
2282 12 : Real h2 = h_from_p_T(_pressure_max, _temperature_min);
2283 12 : Real h3 = h_from_p_T(_pressure_min, _temperature_max);
2284 12 : Real h4 = h_from_p_T(_pressure_max, _temperature_max);
2285 12 : _h_min = std::min({h1, h2, h3, h4});
2286 12 : _h_max = std::max({h1, h2, h3, h4});
2287 : }
2288 : // if csv exists, get max and min values from csv file
2289 38 : else if (_properties.size())
2290 : {
2291 76 : _h_max = *max_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
2292 38 : _h_min = *min_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
2293 : }
2294 0 : else if (_properties_ve.size())
2295 : {
2296 0 : _h_max = *max_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
2297 0 : _h_min = *min_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
2298 : }
2299 : else
2300 0 : mooseError("Need a source to compute the enthalpy grid bounds: either a FP object, or a (p,T) "
2301 : "tabulation file or a (v,e) tabulation file");
2302 :
2303 : // Create h grid for interpolation
2304 : // enthalpy & internal energy use same # grid points
2305 50 : _enthalpy.resize(_num_e);
2306 50 : if (_log_space_h)
2307 : {
2308 : // incrementing the exponent linearly will yield a log-spaced grid after taking the value to
2309 : // the power of 10
2310 0 : if (_h_min < 0)
2311 0 : mooseError("Logarithmic grid in specific energy can only be used with a positive enthalpy. "
2312 0 : "Current minimum: " +
2313 : std::to_string(_h_min));
2314 0 : Real dh = (std::log10(_h_max) - std::log10(_h_min)) / ((Real)_num_e - 1);
2315 0 : Real log_h_min = std::log10(_h_min);
2316 0 : for (const auto j : make_range(_num_e))
2317 0 : _enthalpy[j] = std::pow(10, log_h_min + j * dh);
2318 : }
2319 : else
2320 : {
2321 50 : Real dh = (_h_max - _h_min) / ((Real)_num_e - 1);
2322 3970 : for (const auto j : make_range(_num_e))
2323 3920 : _enthalpy[j] = _h_min + j * dh;
2324 : }
2325 50 : }
2326 :
2327 : void
2328 0 : TabulatedFluidProperties::missingVEInterpolationError(const std::string & function_name) const
2329 : {
2330 0 : mooseError(function_name +
2331 : ": to call this function you must:\n-add this property to the list to the list of "
2332 : "'interpolated_properties'\n and then either:\n-construct (p, T) from (v, e) "
2333 : "tabulations using the 'construct_pT_from_ve' parameter\n-load (v,e) interpolation "
2334 : "tables using the 'fluid_property_ve_file' parameter");
2335 : }
2336 :
2337 : template void TabulatedFluidProperties::checkInputVariables(Real & pressure,
2338 : Real & temperature) const;
2339 : template void TabulatedFluidProperties::checkInputVariables(ADReal & pressure,
2340 : ADReal & temperature) const;
2341 : template void TabulatedFluidProperties::checkInputVariablesVE(Real & v, Real & e) const;
2342 : template void TabulatedFluidProperties::checkInputVariablesVE(ADReal & v, ADReal & e) const;
|
