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 "TabulatedBicubicFluidProperties.h"
11 : #include "BicubicInterpolation.h"
12 : #include "MooseUtils.h"
13 : #include "Conversion.h"
14 : #include "KDTree.h"
15 :
16 : // C++ includes
17 : #include <fstream>
18 : #include <ctime>
19 :
20 : registerMooseObject("FluidPropertiesApp", TabulatedBicubicFluidProperties);
21 : registerMooseObjectRenamed("FluidPropertiesApp",
22 : TabulatedFluidProperties,
23 : "01/01/2023 00:00",
24 : TabulatedBicubicFluidProperties);
25 :
26 : InputParameters
27 597 : TabulatedBicubicFluidProperties::validParams()
28 : {
29 597 : InputParameters params = TabulatedFluidProperties::validParams();
30 597 : params.addClassDescription(
31 : "Fluid properties using bicubic interpolation on tabulated values provided");
32 597 : return params;
33 0 : }
34 :
35 326 : TabulatedBicubicFluidProperties::TabulatedBicubicFluidProperties(const InputParameters & parameters)
36 326 : : TabulatedFluidProperties(parameters)
37 : {
38 308 : }
39 :
40 : void
41 188 : TabulatedBicubicFluidProperties::constructInterpolation()
42 : {
43 : // Construct bicubic interpolants from tabulated data
44 : std::vector<std::vector<Real>> data_matrix;
45 :
46 188 : if (_create_direct_pT_interpolations)
47 : {
48 164 : _property_ipol.resize(_properties.size());
49 1522 : for (const auto i : index_range(_property_ipol))
50 : {
51 1358 : reshapeData2D(_num_p, _num_T, _properties[i], data_matrix);
52 : _property_ipol[i] =
53 1358 : std::make_unique<BicubicInterpolation>(_pressure, _temperature, data_matrix);
54 : }
55 : }
56 :
57 188 : if (_create_direct_ve_interpolations)
58 : {
59 24 : _property_ve_ipol.resize(_properties_ve.size());
60 :
61 264 : for (const auto i : index_range(_property_ve_ipol))
62 : {
63 240 : reshapeData2D(_num_v, _num_e, _properties_ve[i], data_matrix);
64 : _property_ve_ipol[i] =
65 240 : std::make_unique<BicubicInterpolation>(_specific_volume, _internal_energy, data_matrix);
66 : }
67 : }
68 :
69 188 : bool conversion_succeeded = true;
70 : unsigned int fail_counter_ve = 0;
71 : unsigned int fail_counter_vh = 0;
72 : // Create interpolations of (p,T) from (v,e)
73 188 : if (_construct_pT_from_ve)
74 : {
75 : // Grids in specific volume and internal energy can be either linear or logarithmic
76 : // NOTE: this could have been called already when generating tabulated data
77 50 : createVEGridVectors();
78 :
79 : // initialize vectors for interpolation
80 50 : std::vector<std::vector<Real>> p_from_v_e(_num_v);
81 50 : std::vector<std::vector<Real>> T_from_v_e(_num_v);
82 :
83 50 : unsigned int num_p_nans_ve = 0, num_T_nans_ve = 0, num_p_out_bounds_ve = 0,
84 50 : num_T_out_bounds_ve = 0;
85 :
86 3970 : for (unsigned int i = 0; i < _num_v; ++i)
87 : {
88 3920 : Real p_guess = _p_initial_guess;
89 3920 : Real T_guess = _T_initial_guess;
90 3920 : p_from_v_e[i].resize(_num_e);
91 3920 : T_from_v_e[i].resize(_num_e);
92 385120 : for (unsigned int j = 0; j < _num_e; ++j)
93 : {
94 : Real p_ve, T_ve;
95 : // Using an input fluid property instead of the tabulation will get more exact inversions
96 381200 : if (_fp)
97 120000 : _fp->p_T_from_v_e(_specific_volume[i],
98 120000 : _internal_energy[j],
99 120000 : _p_initial_guess,
100 120000 : _T_initial_guess,
101 : p_ve,
102 : T_ve,
103 : conversion_succeeded);
104 : else
105 : {
106 : // The inversion may step outside of the domain of definition of the interpolations,
107 : // which are restricted to the range of the input CSV file
108 261200 : const auto old_error_behavior = _OOBBehavior;
109 261200 : _OOBBehavior = SetToClosestBound;
110 261200 : p_T_from_v_e(_specific_volume[i],
111 261200 : _internal_energy[j],
112 : p_guess,
113 : T_guess,
114 : p_ve,
115 : T_ve,
116 : conversion_succeeded);
117 261200 : _OOBBehavior = old_error_behavior;
118 : // track number of times convergence failed
119 261200 : if (!conversion_succeeded)
120 86110 : ++fail_counter_ve;
121 261200 : }
122 :
123 : // check for NaNs in Newton Method
124 381200 : checkNaNs(_pressure_min,
125 : _pressure_max,
126 : _temperature_min,
127 : _temperature_max,
128 : i,
129 : p_ve,
130 : T_ve,
131 : num_p_nans_ve,
132 : num_T_nans_ve);
133 :
134 : // replace out of bounds pressure values with pmax or pmin
135 381200 : checkOutofBounds(_pressure_min, _pressure_max, p_ve, num_p_out_bounds_ve);
136 : // replace out of bounds temperature values with Tmax or Tmin
137 381200 : checkOutofBounds(_temperature_min, _temperature_max, T_ve, num_T_out_bounds_ve);
138 :
139 381200 : p_from_v_e[i][j] = p_ve;
140 381200 : T_from_v_e[i][j] = T_ve;
141 :
142 381200 : p_guess = p_ve;
143 : T_guess = T_ve;
144 : }
145 : }
146 : // output warning if nans or values out of bounds
147 50 : outputWarnings(num_p_nans_ve,
148 : num_T_nans_ve,
149 : num_p_out_bounds_ve,
150 : num_T_out_bounds_ve,
151 : fail_counter_ve,
152 50 : _num_e * _num_v,
153 : "(v,e)");
154 :
155 : // the bicubic interpolation object are init'ed now
156 : _p_from_v_e_ipol =
157 50 : std::make_unique<BicubicInterpolation>(_specific_volume, _internal_energy, p_from_v_e);
158 : _T_from_v_e_ipol =
159 50 : std::make_unique<BicubicInterpolation>(_specific_volume, _internal_energy, T_from_v_e);
160 50 : }
161 :
162 : // Create interpolations of (p,T) from (v,h)
163 188 : if (_construct_pT_from_vh)
164 : {
165 : // Grids in specific volume and enthalpy can be either linear or logarithmic
166 : // NOTE: the specific volume grid could have been created when generating tabulated data
167 50 : createVHGridVectors();
168 :
169 : // initialize vectors for interpolation
170 50 : std::vector<std::vector<Real>> p_from_v_h(_num_v);
171 50 : std::vector<std::vector<Real>> T_from_v_h(_num_v);
172 :
173 50 : unsigned int num_p_nans_vh = 0, num_T_nans_vh = 0, num_p_out_bounds_vh = 0,
174 50 : num_T_out_bounds_vh = 0;
175 :
176 3970 : for (unsigned int i = 0; i < _num_v; ++i)
177 : {
178 3920 : Real p_guess = _p_initial_guess;
179 3920 : Real T_guess = _T_initial_guess;
180 3920 : p_from_v_h[i].resize(_num_e);
181 3920 : T_from_v_h[i].resize(_num_e);
182 385120 : for (unsigned int j = 0; j < _num_e; ++j)
183 : {
184 : Real p_vh, T_vh;
185 : // Using an input fluid property instead of the tabulation will get more exact inversions
186 381200 : if (_fp)
187 120000 : _fp->p_T_from_v_h(_specific_volume[i],
188 120000 : _enthalpy[j],
189 120000 : _p_initial_guess,
190 120000 : _T_initial_guess,
191 : p_vh,
192 : T_vh,
193 : conversion_succeeded);
194 : else
195 : {
196 : // The inversion may step outside of the domain of definition of the interpolations,
197 : // which are restricted to the range of the input CSV file
198 261200 : const auto old_error_behavior = _OOBBehavior;
199 261200 : _OOBBehavior = SetToClosestBound;
200 261200 : p_T_from_v_h(_specific_volume[i],
201 261200 : _enthalpy[j],
202 : p_guess,
203 : T_guess,
204 : p_vh,
205 : T_vh,
206 : conversion_succeeded);
207 261200 : _OOBBehavior = old_error_behavior;
208 : // track number of times convergence failed
209 261200 : if (!conversion_succeeded)
210 88292 : ++fail_counter_vh;
211 261200 : }
212 :
213 : // check for NaNs in Newton Method
214 381200 : checkNaNs(_pressure_min,
215 : _pressure_max,
216 : _temperature_min,
217 : _temperature_max,
218 : i,
219 : p_vh,
220 : T_vh,
221 : num_p_nans_vh,
222 : num_T_nans_vh);
223 :
224 : // replace out of bounds pressure values with pmax or pmin
225 381200 : checkOutofBounds(_pressure_min, _pressure_max, p_vh, num_p_out_bounds_vh);
226 : // replace out of bounds temperature values with Tmax or Tmin
227 381200 : checkOutofBounds(_temperature_min, _temperature_max, T_vh, num_T_out_bounds_vh);
228 :
229 381200 : p_from_v_h[i][j] = p_vh;
230 381200 : T_from_v_h[i][j] = T_vh;
231 :
232 381200 : p_guess = p_vh;
233 : T_guess = T_vh;
234 : }
235 : }
236 : // output warnings if nans our values out of bounds
237 50 : outputWarnings(num_p_nans_vh,
238 : num_T_nans_vh,
239 : num_p_out_bounds_vh,
240 : num_T_out_bounds_vh,
241 : fail_counter_vh,
242 50 : _num_e * _num_v,
243 : "(v,h)");
244 :
245 : // the bicubic interpolation object are init'ed now
246 : _p_from_v_h_ipol =
247 50 : std::make_unique<BicubicInterpolation>(_specific_volume, _enthalpy, p_from_v_h);
248 : _T_from_v_h_ipol =
249 50 : std::make_unique<BicubicInterpolation>(_specific_volume, _enthalpy, T_from_v_h);
250 50 : }
251 188 : }
252 :
253 : void
254 1598 : TabulatedBicubicFluidProperties::reshapeData2D(unsigned int nrow,
255 : unsigned int ncol,
256 : const std::vector<Real> & vec,
257 : std::vector<std::vector<Real>> & mat)
258 : {
259 1598 : if (!vec.empty())
260 : {
261 1598 : mat.resize(nrow);
262 139326 : for (unsigned int i = 0; i < nrow; ++i)
263 137728 : mat[i].resize(ncol);
264 :
265 139326 : for (unsigned int i = 0; i < nrow; ++i)
266 13700416 : for (unsigned int j = 0; j < ncol; ++j)
267 13562688 : mat[i][j] = vec[i * ncol + j];
268 : }
269 1598 : }
270 :
271 : void
272 762400 : TabulatedBicubicFluidProperties::checkNaNs(Real min_1,
273 : Real max_1,
274 : Real min_2,
275 : Real max_2,
276 : unsigned int i,
277 : Real & variable_1,
278 : Real & variable_2,
279 : unsigned int & num_nans_1,
280 : unsigned int & num_nans_2)
281 : {
282 : // replace nan values with pmax or pmin
283 762400 : if (std::isnan(variable_1))
284 : {
285 96 : if (_specific_volume[i] > ((_v_min + _v_max) / 2))
286 96 : variable_1 = min_1;
287 0 : else if (_specific_volume[i] < ((_v_min + _v_max) / 2))
288 0 : variable_1 = max_1;
289 96 : num_nans_1++;
290 : }
291 : // replace nan values with Tmax or Tmin
292 762400 : if (std::isnan(variable_2))
293 : {
294 96 : if (_specific_volume[i] > ((_v_min + _v_max) / 2))
295 96 : variable_2 = max_2;
296 0 : else if (_specific_volume[i] < ((_v_min + _v_max) / 2))
297 0 : variable_2 = min_2;
298 96 : num_nans_2++;
299 : }
300 762400 : }
301 :
302 : void
303 1524800 : TabulatedBicubicFluidProperties::checkOutofBounds(Real min,
304 : Real max,
305 : Real & variable,
306 : unsigned int & num_out_bounds)
307 : {
308 1524800 : if (variable < min)
309 : {
310 221446 : variable = min;
311 221446 : num_out_bounds++;
312 : }
313 1303354 : else if (variable > max)
314 : {
315 73930 : variable = max;
316 73930 : num_out_bounds++;
317 : }
318 1524800 : }
319 :
320 : void
321 100 : TabulatedBicubicFluidProperties::outputWarnings(unsigned int num_nans_p,
322 : unsigned int num_nans_T,
323 : unsigned int num_out_bounds_p,
324 : unsigned int num_out_bounds_T,
325 : unsigned int convergence_failures,
326 : unsigned int number_points,
327 : std::string variable_set)
328 : {
329 : // make string variables before mooseWarning
330 : std::string while_creating =
331 100 : "While creating (p,T) from " + variable_set + " interpolation tables,\n";
332 100 : std::string warning_message = while_creating;
333 300 : std::string converge_fails = "Inversion to (p,T) from " + variable_set + " failed " +
334 100 : std::to_string(convergence_failures) + " times\n";
335 300 : std::string p_nans = "- " + std::to_string(num_nans_p) + " nans generated out of " +
336 100 : std::to_string(number_points) + " points for pressure\n";
337 300 : std::string T_nans = "- " + std::to_string(num_nans_T) + " nans generated out of " +
338 100 : std::to_string(number_points) + " points for temperature\n";
339 300 : std::string p_oob = "- " + std::to_string(num_out_bounds_p) + " of " +
340 100 : std::to_string(number_points) + " pressure values were out of bounds\n";
341 300 : std::string T_oob = "- " + std::to_string(num_out_bounds_T) + " of " +
342 100 : std::to_string(number_points) + " temperature values were out of bounds\n";
343 : std::string outcome = "The pressure and temperature values were replaced with their respective "
344 100 : "min and max values.\n";
345 :
346 : // bounds are different depending on how the object is used
347 176 : std::string source = (_fp ? "from input parameters" : "from tabulation file");
348 :
349 : // if any of these do not exist, do not want to print them
350 100 : if (convergence_failures)
351 : warning_message += converge_fails;
352 100 : if (num_nans_p)
353 : warning_message += p_nans;
354 100 : if (num_nans_T)
355 : warning_message += T_nans;
356 100 : if (num_out_bounds_p)
357 : {
358 : warning_message += p_oob;
359 400 : warning_message += ("Pressure bounds " + source + ": [" + std::to_string(_pressure_min) + ", " +
360 200 : std::to_string(_pressure_max) + "]\n");
361 : }
362 100 : if (num_out_bounds_T)
363 : {
364 : warning_message += T_oob;
365 300 : warning_message += ("Temperature bounds " + source + ": [" + std::to_string(_temperature_min) +
366 300 : ", " + std::to_string(_temperature_max) + "]\n");
367 : }
368 : // print warning
369 100 : if (num_nans_p || num_nans_T || num_out_bounds_p || num_out_bounds_T || convergence_failures)
370 200 : mooseWarning(warning_message + outcome);
371 100 : }
|
