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