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 1 : if (_fp)
1457 0 : return _fp->T_from_p_h(pressure, enthalpy);
1458 : else
1459 : {
1460 : auto lambda =
1461 1 : [&](ADReal pressure, ADReal current_T, ADReal & new_h, ADReal & dh_dp, ADReal & dh_dT)
1462 : {
1463 1 : h_from_p_T(pressure.value(), current_T.value(), new_h.value(), dh_dp.value(), dh_dT.value());
1464 : // Reconstruct derivatives
1465 : new_h.derivatives() =
1466 1 : dh_dp.value() * pressure.derivatives() + dh_dT.value() * current_T.derivatives();
1467 1 : };
1468 : ADReal T =
1469 2 : FluidPropertiesUtils::NewtonSolve(
1470 2 : pressure, enthalpy, _T_initial_guess, _tolerance, lambda, name() + "::T_from_p_h")
1471 1 : .first;
1472 : // check for nans
1473 1 : if (std::isnan(T))
1474 0 : mooseError("Conversion from enthalpy (h = ",
1475 : enthalpy,
1476 : ") and pressure (p = ",
1477 : pressure,
1478 : ") to temperature failed to converge.");
1479 1 : return T;
1480 : }
1481 : }
1482 :
1483 : Real
1484 0 : TabulatedFluidProperties::s_from_h_p(Real enthalpy, Real pressure) const
1485 : {
1486 0 : Real T = T_from_p_h(pressure, enthalpy);
1487 0 : return s_from_p_T(pressure, T);
1488 : }
1489 :
1490 : void
1491 0 : TabulatedFluidProperties::s_from_h_p(
1492 : Real h, Real pressure, Real & s, Real & ds_dh, Real & ds_dp) const
1493 : {
1494 :
1495 0 : if (_fp)
1496 0 : _fp->s_from_h_p(h, pressure, s, ds_dh, ds_dp);
1497 : else
1498 0 : mooseError("fp", "s_from_h_p derivatives not implemented.");
1499 0 : }
1500 :
1501 : [[noreturn]] void
1502 0 : TabulatedFluidProperties::FluidPropertiesForwardError(const std::string & desired_routine) const
1503 : {
1504 0 : mooseError("TabulatedFluidProperties can only call the function '" + desired_routine +
1505 : "' when the 'fp' parameter is provided. It is currently not implemented using "
1506 : "tabulations, and this property is simply forwarded to the FluidProperties specified "
1507 : "in the 'fp' parameter");
1508 : }
1509 :
1510 : void
1511 24 : TabulatedFluidProperties::writeTabulatedData(std::string file_name)
1512 : {
1513 48 : file_name = file_name.empty() ? "fluid_properties_" + name() + "_out.csv" : file_name;
1514 24 : if (processor_id() == 0)
1515 : {
1516 : {
1517 16 : MooseUtils::checkFileWriteable(file_name);
1518 :
1519 16 : std::ofstream file_out(file_name.c_str());
1520 :
1521 : // Write out date and fluid type
1522 16 : time_t now = std::time(&now);
1523 16 : if (_fp)
1524 16 : file_out << "# " << _fp->fluidName()
1525 16 : << " properties created by TabulatedFluidProperties on " << ctime(&now) << "\n";
1526 : else
1527 8 : file_out << "# tabulated properties created by TabulatedFluidProperties on " << ctime(&now)
1528 16 : << "\n";
1529 :
1530 : // Write out column names
1531 16 : file_out << "pressure, temperature";
1532 160 : for (std::size_t i = 0; i < _interpolated_properties.size(); ++i)
1533 144 : file_out << ", " << _interpolated_properties[i];
1534 16 : file_out << "\n";
1535 :
1536 : // Write out the fluid property data
1537 1616 : for (unsigned int p = 0; p < _num_p; ++p)
1538 161600 : for (unsigned int t = 0; t < _num_T; ++t)
1539 : {
1540 160000 : file_out << _pressure[p] << ", " << _temperature[t];
1541 1600000 : for (std::size_t i = 0; i < _properties.size(); ++i)
1542 1440000 : file_out << ", " << _properties[i][p * _num_T + t];
1543 160000 : file_out << "\n";
1544 : }
1545 :
1546 : file_out << std::flush;
1547 16 : file_out.close();
1548 16 : }
1549 :
1550 : // Write out the (v,e) to (p,T) conversions
1551 16 : if (_construct_pT_from_ve)
1552 : {
1553 8 : const auto file_name_ve = (_file_name_ve_out == "")
1554 16 : ? std::regex_replace(file_name, std::regex("\\.csv"), "_ve.csv")
1555 8 : : _file_name_ve_out;
1556 8 : MooseUtils::checkFileWriteable(file_name_ve);
1557 8 : std::ofstream file_out(file_name_ve.c_str());
1558 :
1559 : // Write out column names
1560 8 : file_out << "specific_volume, internal_energy, pressure, temperature";
1561 80 : for (const auto i : index_range(_properties))
1562 72 : if (_interpolated_properties[i] != "internal_energy")
1563 64 : file_out << ", " << _interpolated_properties[i];
1564 8 : file_out << "\n";
1565 :
1566 : // Write out the fluid property data
1567 88 : for (const auto v : make_range(_num_v))
1568 880 : for (const auto e : make_range(_num_e))
1569 : {
1570 800 : const auto v_val = _specific_volume[v];
1571 800 : const auto e_val = _internal_energy[e];
1572 800 : const auto pressure = _p_from_v_e_ipol->sample(v_val, e_val);
1573 800 : const auto temperature = _T_from_v_e_ipol->sample(v_val, e_val);
1574 3200 : file_out << v_val << ", " << e_val << ", " << pressure << ", " << temperature << ", ";
1575 8000 : for (const auto i : index_range(_properties))
1576 : {
1577 : bool add_comma = true;
1578 7200 : if (i == _density_idx)
1579 800 : file_out << 1 / v_val;
1580 6400 : else if (i == _enthalpy_idx)
1581 800 : file_out << h_from_p_T(pressure, temperature);
1582 : // Note that we could use (p,T) routine to generate this instead of (v,e)
1583 : // Or could use the _properties_ve array
1584 5600 : else if (i == _viscosity_idx)
1585 800 : file_out << mu_from_v_e(v_val, e_val);
1586 4800 : else if (i == _k_idx)
1587 800 : file_out << k_from_v_e(v_val, e_val);
1588 4000 : else if (i == _c_idx)
1589 800 : file_out << c_from_v_e(v_val, e_val);
1590 3200 : else if (i == _cv_idx)
1591 800 : file_out << cv_from_v_e(v_val, e_val);
1592 2400 : else if (i == _cp_idx)
1593 800 : file_out << cp_from_v_e(v_val, e_val);
1594 1600 : else if (i == _entropy_idx)
1595 800 : file_out << s_from_v_e(v_val, e_val);
1596 : else
1597 : add_comma = false;
1598 7200 : if (i != _properties.size() - 1 && add_comma)
1599 5600 : file_out << ", ";
1600 : }
1601 :
1602 800 : file_out << "\n";
1603 : }
1604 :
1605 : file_out << std::flush;
1606 8 : file_out.close();
1607 8 : }
1608 : }
1609 24 : }
1610 :
1611 : void
1612 126 : TabulatedFluidProperties::generateTabulatedData()
1613 : {
1614 : mooseAssert(_fp, "We should not try to generate (p,T) tabulated data without a _fp user object");
1615 126 : _pressure.resize(_num_p);
1616 126 : _temperature.resize(_num_T);
1617 :
1618 : // Generate data for all properties entered in input file
1619 126 : _properties.resize(_interpolated_properties_enum.size());
1620 126 : _interpolated_properties.resize(_interpolated_properties_enum.size());
1621 :
1622 1144 : for (std::size_t i = 0; i < _interpolated_properties_enum.size(); ++i)
1623 1018 : _interpolated_properties[i] = _interpolated_properties_enum[i];
1624 :
1625 1144 : for (const auto i : index_range(_properties))
1626 1018 : _properties[i].resize(_num_p * _num_T);
1627 :
1628 : // Temperature is divided equally into _num_T segments
1629 126 : Real delta_T = (_temperature_max - _temperature_min) / static_cast<Real>(_num_T - 1);
1630 :
1631 12632 : for (unsigned int j = 0; j < _num_T; ++j)
1632 12506 : _temperature[j] = _temperature_min + j * delta_T;
1633 :
1634 : // Divide the pressure into _num_p equal segments
1635 126 : Real delta_p = (_pressure_max - _pressure_min) / static_cast<Real>(_num_p - 1);
1636 :
1637 12632 : for (unsigned int i = 0; i < _num_p; ++i)
1638 12506 : _pressure[i] = _pressure_min + i * delta_p;
1639 :
1640 : // Generate the tabulated data at the pressure and temperature points
1641 1144 : for (const auto i : index_range(_properties))
1642 : {
1643 1018 : if (_interpolated_properties[i] == "density")
1644 12632 : for (unsigned int p = 0; p < _num_p; ++p)
1645 1262542 : for (unsigned int t = 0; t < _num_T; ++t)
1646 1250036 : _properties[i][p * _num_T + t] = _fp->rho_from_p_T(_pressure[p], _temperature[t]);
1647 :
1648 1018 : if (_interpolated_properties[i] == "enthalpy")
1649 12632 : for (unsigned int p = 0; p < _num_p; ++p)
1650 1262542 : for (unsigned int t = 0; t < _num_T; ++t)
1651 1250036 : _properties[i][p * _num_T + t] = _fp->h_from_p_T(_pressure[p], _temperature[t]);
1652 :
1653 1018 : if (_interpolated_properties[i] == "internal_energy")
1654 12632 : for (unsigned int p = 0; p < _num_p; ++p)
1655 1262542 : for (unsigned int t = 0; t < _num_T; ++t)
1656 1250036 : _properties[i][p * _num_T + t] = _fp->e_from_p_T(_pressure[p], _temperature[t]);
1657 :
1658 1018 : if (_interpolated_properties[i] == "viscosity")
1659 12632 : for (unsigned int p = 0; p < _num_p; ++p)
1660 1262542 : for (unsigned int t = 0; t < _num_T; ++t)
1661 1250036 : _properties[i][p * _num_T + t] = _fp->mu_from_p_T(_pressure[p], _temperature[t]);
1662 :
1663 1018 : if (_interpolated_properties[i] == "k")
1664 10309 : for (unsigned int p = 0; p < _num_p; ++p)
1665 1030242 : for (unsigned int t = 0; t < _num_T; ++t)
1666 1020036 : _properties[i][p * _num_T + t] = _fp->k_from_p_T(_pressure[p], _temperature[t]);
1667 :
1668 1018 : if (_interpolated_properties[i] == "c")
1669 10302 : for (unsigned int p = 0; p < _num_p; ++p)
1670 1030200 : for (unsigned int t = 0; t < _num_T; ++t)
1671 1020000 : _properties[i][p * _num_T + t] = _fp->c_from_p_T(_pressure[p], _temperature[t]);
1672 :
1673 1018 : if (_interpolated_properties[i] == "cv")
1674 10309 : for (unsigned int p = 0; p < _num_p; ++p)
1675 1030242 : for (unsigned int t = 0; t < _num_T; ++t)
1676 1020036 : _properties[i][p * _num_T + t] = _fp->cv_from_p_T(_pressure[p], _temperature[t]);
1677 :
1678 1018 : if (_interpolated_properties[i] == "cp")
1679 10309 : for (unsigned int p = 0; p < _num_p; ++p)
1680 1030242 : for (unsigned int t = 0; t < _num_T; ++t)
1681 1020036 : _properties[i][p * _num_T + t] = _fp->cp_from_p_T(_pressure[p], _temperature[t]);
1682 :
1683 1018 : if (_interpolated_properties[i] == "entropy")
1684 10309 : for (unsigned int p = 0; p < _num_p; ++p)
1685 1030242 : for (unsigned int t = 0; t < _num_T; ++t)
1686 1020036 : _properties[i][p * _num_T + t] = _fp->s_from_p_T(_pressure[p], _temperature[t]);
1687 : }
1688 126 : }
1689 :
1690 : void
1691 2 : TabulatedFluidProperties::generateVETabulatedData()
1692 : {
1693 : mooseAssert(_fp, "We should not try to generate (v,e) tabulated data without a _fp user object");
1694 2 : _specific_volume.resize(_num_v);
1695 2 : _internal_energy.resize(_num_e);
1696 :
1697 : // Generate data for all properties entered in input file
1698 2 : _properties_ve.resize(_interpolated_properties_enum.size());
1699 2 : _interpolated_properties.resize(_interpolated_properties_enum.size());
1700 :
1701 : // This is filled from the user input, so it does not matter than this operation is performed
1702 : // for both the (p,T) and (v,e) tabulated data generation
1703 22 : for (std::size_t i = 0; i < _interpolated_properties_enum.size(); ++i)
1704 20 : _interpolated_properties[i] = _interpolated_properties_enum[i];
1705 :
1706 22 : for (const auto i : index_range(_properties_ve))
1707 20 : _properties_ve[i].resize(_num_v * _num_e);
1708 :
1709 : // Grids in (v,e) are not read, so we either use user input or rely on (p,T) data
1710 2 : createVEGridVectors();
1711 :
1712 : // Generate the tabulated data at the pressure and temperature points
1713 22 : for (const auto i : index_range(_properties_ve))
1714 : {
1715 20 : if (_interpolated_properties[i] == "density")
1716 202 : for (unsigned int v = 0; v < _num_v; ++v)
1717 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1718 20000 : _properties_ve[i][v * _num_e + e] = 1. / _specific_volume[v];
1719 :
1720 20 : if (_interpolated_properties[i] == "enthalpy")
1721 202 : for (unsigned int v = 0; v < _num_v; ++v)
1722 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1723 20000 : _properties_ve[i][v * _num_e + e] =
1724 20000 : _fp->h_from_v_e(_specific_volume[v], _internal_energy[e]);
1725 :
1726 20 : if (_interpolated_properties[i] == "internal_energy")
1727 0 : for (unsigned int v = 0; v < _num_v; ++v)
1728 0 : for (unsigned int e = 0; e < _num_e; ++e)
1729 0 : _properties_ve[i][v * _num_e + e] = _internal_energy[e];
1730 :
1731 20 : if (_interpolated_properties[i] == "viscosity")
1732 202 : for (unsigned int v = 0; v < _num_v; ++v)
1733 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1734 20000 : _properties_ve[i][v * _num_e + e] =
1735 20000 : _fp->mu_from_v_e(_specific_volume[v], _internal_energy[e]);
1736 :
1737 20 : if (_interpolated_properties[i] == "k")
1738 202 : for (unsigned int v = 0; v < _num_v; ++v)
1739 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1740 20000 : _properties_ve[i][v * _num_e + e] =
1741 20000 : _fp->k_from_v_e(_specific_volume[v], _internal_energy[e]);
1742 :
1743 20 : if (_interpolated_properties[i] == "c")
1744 202 : for (unsigned int v = 0; v < _num_v; ++v)
1745 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1746 20000 : _properties_ve[i][v * _num_e + e] =
1747 20000 : _fp->c_from_v_e(_specific_volume[v], _internal_energy[e]);
1748 :
1749 20 : if (_interpolated_properties[i] == "cv")
1750 202 : for (unsigned int v = 0; v < _num_v; ++v)
1751 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1752 20000 : _properties_ve[i][v * _num_e + e] =
1753 20000 : _fp->cv_from_v_e(_specific_volume[v], _internal_energy[e]);
1754 :
1755 20 : if (_interpolated_properties[i] == "cp")
1756 202 : for (unsigned int v = 0; v < _num_v; ++v)
1757 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1758 20000 : _properties_ve[i][v * _num_e + e] =
1759 20000 : _fp->cp_from_v_e(_specific_volume[v], _internal_energy[e]);
1760 :
1761 20 : if (_interpolated_properties[i] == "entropy")
1762 202 : for (unsigned int v = 0; v < _num_v; ++v)
1763 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1764 20000 : _properties_ve[i][v * _num_e + e] =
1765 20000 : _fp->s_from_v_e(_specific_volume[v], _internal_energy[e]);
1766 :
1767 20 : if (_interpolated_properties[i] == "pressure")
1768 202 : for (unsigned int v = 0; v < _num_v; ++v)
1769 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1770 20000 : _properties_ve[i][v * _num_e + e] =
1771 20000 : _fp->p_from_v_e(_specific_volume[v], _internal_energy[e]);
1772 :
1773 20 : if (_interpolated_properties[i] == "temperature")
1774 202 : for (unsigned int v = 0; v < _num_v; ++v)
1775 20200 : for (unsigned int e = 0; e < _num_e; ++e)
1776 20000 : _properties_ve[i][v * _num_e + e] =
1777 20000 : _fp->T_from_v_e(_specific_volume[v], _internal_energy[e]);
1778 : }
1779 2 : }
1780 :
1781 : template <typename T>
1782 : void
1783 37000049 : TabulatedFluidProperties::checkInputVariables(T & pressure, T & temperature) const
1784 : {
1785 37000049 : if (_OOBBehavior == Ignore)
1786 : return;
1787 36999401 : else if (MooseUtils::absoluteFuzzyGreaterThan(_pressure_min, pressure, libMesh::TOLERANCE) ||
1788 : MooseUtils::absoluteFuzzyGreaterThan(pressure, _pressure_max, libMesh::TOLERANCE))
1789 : {
1790 28249314 : if (_OOBBehavior == Throw)
1791 12 : throw MooseException("Pressure " + Moose::stringify(pressure) +
1792 : " is outside the range of tabulated pressure (" +
1793 2 : Moose::stringify(_pressure_min) + ", " +
1794 2 : Moose::stringify(_pressure_max) + ").");
1795 :
1796 : else
1797 : {
1798 36915376 : pressure = std::max(_pressure_min, std::min(pressure, _pressure_max));
1799 28249312 : if (_OOBBehavior == DeclareInvalid)
1800 146 : flagInvalidSolution("Pressure out of bounds");
1801 28249168 : else if (_OOBBehavior == WarnInvalid)
1802 545 : flagSolutionWarning("Pressure out of bounds");
1803 : }
1804 : }
1805 :
1806 36999399 : if (MooseUtils::absoluteFuzzyGreaterThan(_temperature_min, temperature, libMesh::TOLERANCE) ||
1807 : MooseUtils::absoluteFuzzyGreaterThan(temperature, _temperature_max, libMesh::TOLERANCE))
1808 : {
1809 11054558 : if (_OOBBehavior == Throw)
1810 12 : throw MooseException("Temperature " + Moose::stringify(temperature) +
1811 : " is outside the range of tabulated temperature (" +
1812 2 : Moose::stringify(_temperature_min) + ", " +
1813 2 : Moose::stringify(_temperature_max) + ").");
1814 : else
1815 : {
1816 12574348 : temperature = std::max(T(_temperature_min), std::min(temperature, T(_temperature_max)));
1817 11054556 : if (_OOBBehavior == DeclareInvalid)
1818 146 : flagInvalidSolution("Temperature out of bounds");
1819 11054412 : else if (_OOBBehavior == WarnInvalid)
1820 545 : flagSolutionWarning("Temperature out of bounds");
1821 : }
1822 : }
1823 : }
1824 :
1825 : template <typename T>
1826 : void
1827 5285 : TabulatedFluidProperties::checkInputVariablesVE(T & v, T & e) const
1828 : {
1829 5285 : if (_OOBBehavior == Ignore)
1830 : return;
1831 5285 : else if (e < _e_min || e > _e_max)
1832 : {
1833 67 : if (_OOBBehavior == Throw)
1834 0 : throw MooseException("Specific internal energy " + Moose::stringify(e) +
1835 : " is outside the range of tabulated specific internal energies (" +
1836 0 : Moose::stringify(_e_min) + ", " + Moose::stringify(_e_max) + ").");
1837 : else
1838 : {
1839 80 : e = std::max(_e_min, std::min(e, _e_max));
1840 67 : if (_OOBBehavior == DeclareInvalid)
1841 0 : flagInvalidSolution("Specific internal energy out of bounds");
1842 67 : else if (_OOBBehavior == WarnInvalid)
1843 0 : flagSolutionWarning("Specific internal energy out of bounds");
1844 : }
1845 : }
1846 :
1847 5285 : if (v < _v_min || v > _v_max)
1848 : {
1849 54 : if (_OOBBehavior == Throw)
1850 0 : throw MooseException("Specific volume " + Moose::stringify(v) +
1851 : " is outside the range of tabulated specific volumes (" +
1852 0 : Moose::stringify(_v_min) + ", " + Moose::stringify(_v_max) + ").");
1853 : else
1854 : {
1855 54 : v = std::max(T(_v_min), std::min(v, T(_v_max)));
1856 54 : if (_OOBBehavior == DeclareInvalid)
1857 0 : flagInvalidSolution("Specific volume out of bounds");
1858 54 : else if (_OOBBehavior == WarnInvalid)
1859 0 : flagSolutionWarning("Specific volume out of bounds");
1860 : }
1861 : }
1862 : }
1863 :
1864 : void
1865 326 : TabulatedFluidProperties::checkInitialGuess() const
1866 : {
1867 326 : if (_fp && (_construct_pT_from_ve || _construct_pT_from_vh))
1868 : {
1869 16 : if (_p_initial_guess < _pressure_min || _p_initial_guess > _pressure_max)
1870 28 : mooseWarning("Pressure initial guess for (p,T), (v,e) conversions " +
1871 26 : Moose::stringify(_p_initial_guess) +
1872 : " is outside the range of tabulated "
1873 14 : "pressure (" +
1874 40 : Moose::stringify(_pressure_min) + ", " + Moose::stringify(_pressure_max) + ").");
1875 :
1876 14 : if (_T_initial_guess < _temperature_min || _T_initial_guess > _temperature_max)
1877 28 : mooseWarning("Temperature initial guess for (p,T), (v,e) conversions " +
1878 26 : Moose::stringify(_T_initial_guess) +
1879 : " is outside the range of tabulated "
1880 14 : "temperature (" +
1881 40 : Moose::stringify(_temperature_min) + ", " + Moose::stringify(_temperature_max) +
1882 : ").");
1883 : }
1884 322 : }
1885 :
1886 : void
1887 65 : TabulatedFluidProperties::readFileTabulationData(const bool use_pT)
1888 : {
1889 : std::string file_name;
1890 65 : if (use_pT)
1891 : {
1892 43 : _console << name() + ": Reading tabulated properties from " << _file_name_in << std::endl;
1893 43 : _csv_reader.read();
1894 43 : file_name = _file_name_in;
1895 : }
1896 : else
1897 : {
1898 27 : _console << name() + ": Reading tabulated properties from " << _file_name_ve_in << std::endl;
1899 22 : _csv_reader.setFileName(_file_name_ve_in);
1900 22 : _csv_reader.read();
1901 : file_name = _file_name_ve_in;
1902 : }
1903 :
1904 65 : const std::vector<std::string> & column_names = _csv_reader.getNames();
1905 :
1906 : // These columns form the grid and must be present in the file
1907 : std::vector<std::string> required_columns;
1908 65 : if (use_pT)
1909 129 : required_columns = {"pressure", "temperature"};
1910 : else
1911 66 : required_columns = {"specific_volume", "internal_energy"};
1912 :
1913 : // Check that all required columns are present
1914 194 : for (std::size_t i = 0; i < required_columns.size(); ++i)
1915 : {
1916 130 : if (std::find(column_names.begin(), column_names.end(), required_columns[i]) ==
1917 : column_names.end())
1918 1 : mooseError("No ",
1919 : required_columns[i],
1920 : " data read in ",
1921 : file_name,
1922 : ". A column named ",
1923 : required_columns[i],
1924 : " must be present");
1925 : }
1926 :
1927 : // These columns can be present in the file
1928 : std::vector<std::string> property_columns = {
1929 64 : "density", "enthalpy", "viscosity", "k", "c", "cv", "cp", "entropy"};
1930 64 : if (use_pT)
1931 84 : property_columns.push_back("internal_energy");
1932 : else
1933 : {
1934 22 : property_columns.push_back("pressure");
1935 44 : property_columns.push_back("temperature");
1936 : }
1937 :
1938 : // Check that any property names read from the file are present in the list of possible
1939 : // properties, and if they are, add them to the list of read properties
1940 762 : for (std::size_t i = 0; i < column_names.size(); ++i)
1941 : {
1942 : // Only check properties not in _required_columns
1943 699 : if (std::find(required_columns.begin(), required_columns.end(), column_names[i]) ==
1944 : required_columns.end())
1945 : {
1946 571 : if (std::find(property_columns.begin(), property_columns.end(), column_names[i]) ==
1947 : property_columns.end())
1948 1 : mooseWarning(column_names[i],
1949 : " read in ",
1950 : file_name,
1951 : " tabulation file is not one of the properties that TabulatedFluidProperties "
1952 : "understands. It will be ignored.");
1953 : // We could be reading a (v,e) tabulation after having read a (p,T) tabulation, do not
1954 : // insert twice
1955 570 : else if (std::find(_interpolated_properties.begin(),
1956 : _interpolated_properties.end(),
1957 : column_names[i]) == _interpolated_properties.end())
1958 570 : _interpolated_properties.push_back(column_names[i]);
1959 : }
1960 : }
1961 :
1962 : std::map<std::string, unsigned int> data_index;
1963 758 : for (std::size_t i = 0; i < column_names.size(); ++i)
1964 : {
1965 695 : auto it = std::find(column_names.begin(), column_names.end(), column_names[i]);
1966 695 : data_index[column_names[i]] = std::distance(column_names.begin(), it);
1967 : }
1968 :
1969 63 : const std::vector<std::vector<Real>> & column_data = _csv_reader.getData();
1970 :
1971 : // Extract the pressure and temperature data vectors
1972 63 : if (use_pT)
1973 : {
1974 41 : _pressure = column_data[data_index.find("pressure")->second];
1975 82 : _temperature = column_data[data_index.find("temperature")->second];
1976 : }
1977 : else
1978 : {
1979 22 : _specific_volume = column_data[data_index.find("specific_volume")->second];
1980 44 : _internal_energy = column_data[data_index.find("internal_energy")->second];
1981 : }
1982 :
1983 63 : if (use_pT)
1984 82 : checkFileTabulationGrids(_pressure, _temperature, file_name, "pressure", "temperature");
1985 : else
1986 44 : checkFileTabulationGrids(_specific_volume,
1987 22 : _internal_energy,
1988 : file_name,
1989 : "specific volume",
1990 : "specific internal energy");
1991 :
1992 60 : if (use_pT)
1993 : {
1994 38 : _num_p = _pressure.size();
1995 38 : _num_T = _temperature.size();
1996 :
1997 : // Minimum and maximum pressure and temperature. Note that _pressure and
1998 : // _temperature are sorted
1999 38 : _pressure_min = _pressure.front();
2000 38 : _pressure_max = _pressure.back();
2001 38 : _temperature_min = _temperature.front();
2002 38 : _temperature_max = _temperature.back();
2003 :
2004 : // Extract the fluid property data from the file
2005 378 : for (std::size_t i = 0; i < _interpolated_properties.size(); ++i)
2006 340 : _properties.push_back(column_data[data_index.find(_interpolated_properties[i])->second]);
2007 : }
2008 : else
2009 : {
2010 22 : _num_v = _specific_volume.size();
2011 22 : _num_e = _internal_energy.size();
2012 :
2013 : // Minimum and maximum specific internal energy and specific volume
2014 22 : _v_min = _specific_volume.front();
2015 22 : _v_max = _specific_volume.back();
2016 22 : _e_min = _internal_energy.front();
2017 22 : _e_max = _internal_energy.back();
2018 :
2019 : // We cannot overwrite the tabulated data grid with a grid generated from user-input for the
2020 : // purpose of creating (p,T) to (v,e) interpolations
2021 22 : if (_construct_pT_from_ve)
2022 3 : paramError("construct_pT_from_ve",
2023 : "Reading a (v,e) tabulation and generating (p,T) to (v,e) interpolation tables is "
2024 : "not supported at this time.");
2025 :
2026 : // Make sure we use the tabulation bounds
2027 22 : _e_bounds_specified = true;
2028 22 : _v_bounds_specified = true;
2029 :
2030 : // Extract the fluid property data from the file
2031 22 : _properties_ve.reserve(_interpolated_properties.size());
2032 242 : for (std::size_t i = 0; i < _interpolated_properties.size(); ++i)
2033 220 : _properties_ve.push_back(column_data[data_index.find(_interpolated_properties[i])->second]);
2034 : }
2035 259 : }
2036 :
2037 : void
2038 63 : TabulatedFluidProperties::checkFileTabulationGrids(std::vector<Real> & v1,
2039 : std::vector<Real> & v2,
2040 : const std::string & file_name,
2041 : const std::string & v1_name,
2042 : const std::string & v2_name)
2043 : {
2044 : // NOTE: We kept the comments in terms of pressure & temperature for clarity
2045 : // Pressure (v1) and temperature (v2) data contains duplicates due to the csv format.
2046 : // First, check that pressure (v1) is monotonically increasing
2047 63 : if (!std::is_sorted(v1.begin(), v1.end()))
2048 0 : mooseError("The column data for ", v1_name, " is not monotonically increasing in ", file_name);
2049 :
2050 : // The first pressure (v1) value is repeated for each temperature (v2) value. Counting the
2051 : // number of repeats provides the number of temperature (v2) values
2052 63 : auto num_v2 = std::count(v1.begin(), v1.end(), v1.front());
2053 :
2054 : // Now remove the duplicates in the pressure (v1) vector
2055 63 : auto last_unique = std::unique(v1.begin(), v1.end());
2056 : v1.erase(last_unique, v1.end());
2057 :
2058 : // Check that the number of rows in the csv file is equal to _num_v1 * _num_v2
2059 : // v2 is currently the same size as the column_data (will get trimmed at the end)
2060 63 : if (v2.size() != v1.size() * libMesh::cast_int<unsigned int>(num_v2))
2061 1 : mooseError("The number of rows in ",
2062 : file_name,
2063 : " is not equal to the number of unique ",
2064 : v1_name,
2065 : " values ",
2066 1 : v1.size(),
2067 : " multiplied by the number of unique ",
2068 : v2_name,
2069 : " values ",
2070 : num_v2);
2071 :
2072 : // Need to make sure that the temperature (v2) values are provided in ascending order
2073 62 : std::vector<Real> base_v2(v2.begin(), v2.begin() + num_v2);
2074 62 : if (!std::is_sorted(base_v2.begin(), base_v2.end()))
2075 1 : mooseError("The column data for ", v2_name, " is not monotonically increasing in ", file_name);
2076 :
2077 : // Need to make sure that the temperature (v2) are repeated for each pressure (v1) grid point
2078 61 : auto it_v2 = v2.begin() + num_v2;
2079 3832 : for (const auto i : make_range(v1.size() - 1))
2080 : {
2081 3773 : std::vector<Real> repeated_v2(it_v2, it_v2 + num_v2);
2082 3772 : if (repeated_v2 != base_v2)
2083 1 : mooseError(v2_name,
2084 : " values for ",
2085 : v1_name,
2086 : " ",
2087 1 : v1[i + 1],
2088 : " are not identical to values for ",
2089 : v1[0]);
2090 :
2091 : std::advance(it_v2, num_v2);
2092 3772 : }
2093 :
2094 : // At this point, all temperature (v2) data has been provided in ascending order
2095 : // identically for each pressure (v1) value, so we can just keep the first range
2096 : v2.erase(v2.begin() + num_v2, v2.end());
2097 62 : }
2098 :
2099 : void
2100 188 : TabulatedFluidProperties::computePropertyIndicesInInterpolationVectors()
2101 : {
2102 : // At this point, all properties read or generated are able to be used by
2103 : // TabulatedFluidProperties. Now set flags and indexes for each property in
2104 : //_interpolated_properties to use in property calculations
2105 1786 : for (std::size_t i = 0; i < _interpolated_properties.size(); ++i)
2106 : {
2107 1598 : if (_interpolated_properties[i] == "density")
2108 : {
2109 188 : _interpolate_density = true;
2110 188 : _density_idx = i;
2111 : }
2112 1410 : else if (_interpolated_properties[i] == "enthalpy")
2113 : {
2114 188 : _interpolate_enthalpy = true;
2115 188 : _enthalpy_idx = i;
2116 : }
2117 1222 : else if (_interpolated_properties[i] == "internal_energy")
2118 : {
2119 164 : _interpolate_internal_energy = true;
2120 164 : _internal_energy_idx = i;
2121 : }
2122 1058 : else if (_interpolated_properties[i] == "viscosity")
2123 : {
2124 188 : _interpolate_viscosity = true;
2125 188 : _viscosity_idx = i;
2126 : }
2127 870 : else if (_interpolated_properties[i] == "k")
2128 : {
2129 165 : _interpolate_k = true;
2130 165 : _k_idx = i;
2131 : }
2132 705 : else if (_interpolated_properties[i] == "c")
2133 : {
2134 162 : _interpolate_c = true;
2135 162 : _c_idx = i;
2136 : }
2137 543 : else if (_interpolated_properties[i] == "cp")
2138 : {
2139 165 : _interpolate_cp = true;
2140 165 : _cp_idx = i;
2141 : }
2142 378 : else if (_interpolated_properties[i] == "cv")
2143 : {
2144 165 : _interpolate_cv = true;
2145 165 : _cv_idx = i;
2146 : }
2147 213 : else if (_interpolated_properties[i] == "entropy")
2148 : {
2149 165 : _interpolate_entropy = true;
2150 165 : _entropy_idx = i;
2151 : }
2152 48 : else if (_interpolated_properties[i] == "pressure")
2153 : {
2154 24 : _interpolate_pressure = true;
2155 24 : _p_idx = i;
2156 : }
2157 24 : else if (_interpolated_properties[i] == "temperature")
2158 : {
2159 24 : _interpolate_temperature = true;
2160 24 : _T_idx = i;
2161 : }
2162 : else
2163 0 : mooseError("Specified property '" + _interpolated_properties[i] +
2164 : "' is present in the tabulation but is not currently leveraged by the code in the "
2165 : "TabulatedFluidProperties. If it is spelled correctly, then please contact a "
2166 : "MOOSE or fluid properties module developer.");
2167 : }
2168 188 : }
2169 :
2170 : void
2171 102 : TabulatedFluidProperties::createVGridVector()
2172 : {
2173 : mooseAssert(_file_name_ve_in.empty(), "We should be reading the specific volume grid from file");
2174 102 : if (!_v_bounds_specified)
2175 : {
2176 52 : if (_fp)
2177 : {
2178 : // extreme values of specific volume for the grid bounds
2179 14 : Real v1 = v_from_p_T(_pressure_min, _temperature_min);
2180 14 : Real v2 = v_from_p_T(_pressure_max, _temperature_min);
2181 14 : Real v3 = v_from_p_T(_pressure_min, _temperature_max);
2182 14 : Real v4 = v_from_p_T(_pressure_max, _temperature_max);
2183 14 : _v_min = std::min({v1, v2, v3, v4});
2184 14 : _v_max = std::max({v1, v2, v3, v4});
2185 : }
2186 : // if csv exists, get max and min values from csv file
2187 : else
2188 : {
2189 : Real rho_max =
2190 76 : *std::max_element(_properties[_density_idx].begin(), _properties[_density_idx].end());
2191 : Real rho_min =
2192 38 : *std::min_element(_properties[_density_idx].begin(), _properties[_density_idx].end());
2193 38 : _v_max = 1 / rho_min;
2194 38 : _v_min = 1 / rho_max;
2195 : }
2196 : // Prevent changing the bounds of the grid
2197 52 : _v_bounds_specified = true;
2198 : }
2199 :
2200 : // Create v grid for interpolation
2201 102 : _specific_volume.resize(_num_v);
2202 102 : if (_log_space_v)
2203 : {
2204 : // incrementing the exponent linearly will yield a log-spaced grid after taking the value to
2205 : // the power of 10
2206 24 : Real dv = (std::log10(_v_max) - std::log10(_v_min)) / ((Real)_num_v - 1);
2207 24 : Real log_v_min = std::log10(_v_min);
2208 2424 : for (unsigned int j = 0; j < _num_v; ++j)
2209 2400 : _specific_volume[j] = std::pow(10, log_v_min + j * dv);
2210 : }
2211 : else
2212 : {
2213 78 : Real dv = (_v_max - _v_min) / ((Real)_num_v - 1);
2214 5718 : for (unsigned int j = 0; j < _num_v; ++j)
2215 5640 : _specific_volume[j] = _v_min + j * dv;
2216 : }
2217 102 : }
2218 :
2219 : void
2220 52 : TabulatedFluidProperties::createVEGridVectors()
2221 : {
2222 52 : createVGridVector();
2223 52 : if (!_e_bounds_specified)
2224 : {
2225 52 : if (_fp)
2226 : {
2227 : // extreme values of internal energy for the grid bounds
2228 14 : Real e1 = e_from_p_T(_pressure_min, _temperature_min);
2229 14 : Real e2 = e_from_p_T(_pressure_max, _temperature_min);
2230 14 : Real e3 = e_from_p_T(_pressure_min, _temperature_max);
2231 14 : Real e4 = e_from_p_T(_pressure_max, _temperature_max);
2232 14 : _e_min = std::min({e1, e2, e3, e4});
2233 14 : _e_max = std::max({e1, e2, e3, e4});
2234 : }
2235 : // if csv exists, get max and min values from csv file
2236 : else
2237 : {
2238 38 : _e_min = *std::min_element(_properties[_internal_energy_idx].begin(),
2239 38 : _properties[_internal_energy_idx].end());
2240 38 : _e_max = *std::max_element(_properties[_internal_energy_idx].begin(),
2241 : _properties[_internal_energy_idx].end());
2242 : }
2243 : }
2244 :
2245 : // Create e grid for interpolation
2246 52 : _internal_energy.resize(_num_e);
2247 52 : if (_log_space_e)
2248 : {
2249 : // incrementing the exponent linearly will yield a log-spaced grid after taking the value to
2250 : // the power of 10
2251 12 : if (_e_min < 0)
2252 0 : mooseError("Logarithmic grid in specific energy can only be used with a positive specific "
2253 0 : "energy. Current minimum: " +
2254 : std::to_string(_e_min));
2255 12 : Real de = (std::log10(_e_max) - std::log10(_e_min)) / ((Real)_num_e - 1);
2256 12 : Real log_e_min = std::log10(_e_min);
2257 1212 : for (const auto j : make_range(_num_e))
2258 1200 : _internal_energy[j] = std::pow(10, log_e_min + j * de);
2259 : }
2260 : else
2261 : {
2262 40 : Real de = (_e_max - _e_min) / ((Real)_num_e - 1);
2263 2960 : for (const auto j : make_range(_num_e))
2264 2920 : _internal_energy[j] = _e_min + j * de;
2265 : }
2266 52 : }
2267 :
2268 : void
2269 50 : TabulatedFluidProperties::createVHGridVectors()
2270 : {
2271 50 : if (_file_name_ve_in.empty())
2272 50 : createVGridVector();
2273 50 : if (_fp)
2274 : {
2275 : // extreme values of enthalpy for the grid bounds
2276 12 : Real h1 = h_from_p_T(_pressure_min, _temperature_min);
2277 12 : Real h2 = h_from_p_T(_pressure_max, _temperature_min);
2278 12 : Real h3 = h_from_p_T(_pressure_min, _temperature_max);
2279 12 : Real h4 = h_from_p_T(_pressure_max, _temperature_max);
2280 12 : _h_min = std::min({h1, h2, h3, h4});
2281 12 : _h_max = std::max({h1, h2, h3, h4});
2282 : }
2283 : // if csv exists, get max and min values from csv file
2284 38 : else if (_properties.size())
2285 : {
2286 76 : _h_max = *max_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
2287 38 : _h_min = *min_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
2288 : }
2289 0 : else if (_properties_ve.size())
2290 : {
2291 0 : _h_max = *max_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
2292 0 : _h_min = *min_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
2293 : }
2294 : else
2295 0 : mooseError("Need a source to compute the enthalpy grid bounds: either a FP object, or a (p,T) "
2296 : "tabulation file or a (v,e) tabulation file");
2297 :
2298 : // Create h grid for interpolation
2299 : // enthalpy & internal energy use same # grid points
2300 50 : _enthalpy.resize(_num_e);
2301 50 : if (_log_space_h)
2302 : {
2303 : // incrementing the exponent linearly will yield a log-spaced grid after taking the value to
2304 : // the power of 10
2305 0 : if (_h_min < 0)
2306 0 : mooseError("Logarithmic grid in specific energy can only be used with a positive enthalpy. "
2307 0 : "Current minimum: " +
2308 : std::to_string(_h_min));
2309 0 : Real dh = (std::log10(_h_max) - std::log10(_h_min)) / ((Real)_num_e - 1);
2310 0 : Real log_h_min = std::log10(_h_min);
2311 0 : for (const auto j : make_range(_num_e))
2312 0 : _enthalpy[j] = std::pow(10, log_h_min + j * dh);
2313 : }
2314 : else
2315 : {
2316 50 : Real dh = (_h_max - _h_min) / ((Real)_num_e - 1);
2317 3970 : for (const auto j : make_range(_num_e))
2318 3920 : _enthalpy[j] = _h_min + j * dh;
2319 : }
2320 50 : }
2321 :
2322 : void
2323 0 : TabulatedFluidProperties::missingVEInterpolationError(const std::string & function_name) const
2324 : {
2325 0 : mooseError(function_name +
2326 : ": to call this function you must:\n-add this property to the list to the list of "
2327 : "'interpolated_properties'\n and then either:\n-construct (p, T) from (v, e) "
2328 : "tabulations using the 'construct_pT_from_ve' parameter\n-load (v,e) interpolation "
2329 : "tables using the 'fluid_property_ve_file' parameter");
2330 : }
2331 :
2332 : template void TabulatedFluidProperties::checkInputVariables(Real & pressure,
2333 : Real & temperature) const;
2334 : template void TabulatedFluidProperties::checkInputVariables(ADReal & pressure,
2335 : ADReal & temperature) const;
2336 : template void TabulatedFluidProperties::checkInputVariablesVE(Real & v, Real & e) const;
2337 : template void TabulatedFluidProperties::checkInputVariablesVE(ADReal & v, ADReal & e) const;
|
