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 "ViewFactorBase.h"
11 : #include "libmesh/quadrature.h"
12 :
13 : #include <limits>
14 :
15 : InputParameters
16 372 : ViewFactorBase::validParams()
17 : {
18 372 : InputParameters params = SideUserObject::validParams();
19 744 : params.addParam<Real>("view_factor_tol",
20 744 : std::numeric_limits<Real>::max(),
21 : "Tolerance for checking view factors. Default is to allow everything.");
22 744 : params.addParam<bool>(
23 744 : "print_view_factor_info", false, "Flag to print information about computed view factors.");
24 744 : params.addParam<bool>("normalize_view_factor",
25 744 : true,
26 : "Determines if view factors are normalized to sum to one (consistent with "
27 : "their definition).");
28 :
29 744 : MooseEnum normalization_method("lagrange_multiplier inverse_row_sum", "lagrange_multiplier");
30 744 : normalization_method.addDocumentation("lagrange_multiplier", "Uses Lagrange multiplier method");
31 744 : normalization_method.addDocumentation("inverse_row_sum",
32 : "Divides each view factor by its row sum");
33 744 : params.addParam<MooseEnum>("normalization_method",
34 : normalization_method,
35 : "Normalization method to use if 'normalize_view_factor' is 'true'.");
36 :
37 372 : params.addClassDescription(
38 : "A base class for automatic computation of view factors between sidesets.");
39 372 : return params;
40 372 : }
41 :
42 201 : ViewFactorBase::ViewFactorBase(const InputParameters & parameters)
43 : : SideUserObject(parameters),
44 402 : _n_sides(boundaryIDs().size()),
45 201 : _areas(_n_sides),
46 402 : _view_factor_tol(getParam<Real>("view_factor_tol")),
47 402 : _normalize_view_factor(getParam<bool>("normalize_view_factor")),
48 201 : _normalization_method(
49 201 : getParam<MooseEnum>("normalization_method").getEnum<NormalizationMethod>()),
50 603 : _print_view_factor_info(getParam<bool>("print_view_factor_info"))
51 : {
52 : // sizing the view factor array
53 201 : _view_factors.resize(_n_sides);
54 1497 : for (auto & v : _view_factors)
55 1296 : v.resize(_n_sides);
56 :
57 : // set up the map from the side id to the local index & side name to local index
58 402 : std::vector<BoundaryName> boundary_names = getParam<std::vector<BoundaryName>>("boundary");
59 1497 : for (unsigned int j = 0; j < boundary_names.size(); ++j)
60 1296 : _side_name_index[boundary_names[j]] = j;
61 201 : }
62 :
63 : unsigned int
64 7764 : ViewFactorBase::getSideNameIndex(std::string name) const
65 : {
66 : auto it = _side_name_index.find(name);
67 7764 : if (it == _side_name_index.end())
68 0 : mooseError("Boundary ", name, " does not exist.");
69 7764 : return it->second;
70 : }
71 :
72 : Real
73 332 : ViewFactorBase::getViewFactor(BoundaryID from_id, BoundaryID to_id) const
74 : {
75 332 : auto from_name = _mesh.getBoundaryName(from_id);
76 332 : auto to_name = _mesh.getBoundaryName(to_id);
77 :
78 996 : return getViewFactor(from_name, to_name);
79 : }
80 :
81 : Real
82 405611 : ViewFactorBase::getViewFactor(BoundaryName from_name, BoundaryName to_name) const
83 : {
84 405611 : auto from = _side_name_index.find(from_name);
85 405611 : auto to = _side_name_index.find(to_name);
86 405611 : if (from == _side_name_index.end())
87 0 : mooseError("Boundary id ",
88 0 : _mesh.getBoundaryID(from_name),
89 : " with name ",
90 : from_name,
91 : " not listed in boundary parameter.");
92 :
93 405611 : if (to == _side_name_index.end())
94 0 : mooseError("Boundary id ",
95 0 : _mesh.getBoundaryID(to_name),
96 : " with name ",
97 : to_name,
98 : " not listed in boundary parameter.");
99 :
100 405611 : return _view_factors[from->second][to->second];
101 : }
102 :
103 : void
104 164 : ViewFactorBase::initialize()
105 : {
106 : std::fill(_areas.begin(), _areas.end(), 0);
107 164 : }
108 :
109 : void
110 3724 : ViewFactorBase::execute()
111 : {
112 : // compute areas
113 3724 : auto current_boundary_name = _mesh.getBoundaryName(_current_boundary_id);
114 : auto it = _side_name_index.find(current_boundary_name);
115 3724 : if (it == _side_name_index.end())
116 0 : mooseError("Current boundary name: ",
117 : current_boundary_name,
118 : " with id ",
119 0 : _current_boundary_id,
120 : " not in boundary parameter.");
121 :
122 3724 : _areas[it->second] += _current_side_volume;
123 3724 : }
124 :
125 : void
126 142 : ViewFactorBase::finalize()
127 : {
128 142 : gatherSum(_areas);
129 :
130 142 : finalizeViewFactor();
131 142 : checkAndNormalizeViewFactor();
132 142 : }
133 :
134 : void
135 22 : ViewFactorBase::threadJoin(const UserObject & y)
136 : {
137 : const auto & pps = static_cast<const ViewFactorBase &>(y);
138 171 : for (unsigned int i = 0; i < _n_sides; ++i)
139 149 : _areas[i] += pps._areas[i];
140 :
141 22 : threadJoinViewFactor(y);
142 22 : }
143 :
144 : Real
145 14052 : ViewFactorBase::devReciprocity(unsigned int i, unsigned int j) const
146 : {
147 14052 : return _view_factors[i][j] - _areas[j] / _areas[i] * _view_factors[j][i];
148 : }
149 :
150 : Real
151 284 : ViewFactorBase::maxDevReciprocity() const
152 : {
153 : Real v = 0;
154 2184 : for (unsigned int i = 0; i < _n_sides; ++i)
155 7744 : for (unsigned int j = i + 1; j < _n_sides; ++j)
156 : {
157 5844 : Real s = std::abs(devReciprocity(i, j));
158 5844 : if (s > v)
159 : v = s;
160 : }
161 284 : return v;
162 : }
163 :
164 : Real
165 2528 : ViewFactorBase::viewFactorRowSum(unsigned int i) const
166 : {
167 : Real s = 0;
168 20944 : for (unsigned int to = 0; to < _n_sides; ++to)
169 18416 : s += _view_factors[i][to];
170 2528 : return s;
171 : }
172 :
173 : Real
174 284 : ViewFactorBase::maxDevRowSum() const
175 : {
176 : Real v = 0;
177 2184 : for (unsigned int i = 0; i < _n_sides; ++i)
178 : {
179 1900 : Real s = std::abs(1 - viewFactorRowSum(i));
180 1900 : if (s > v)
181 : v = s;
182 : }
183 284 : return v;
184 : }
185 :
186 : void
187 142 : ViewFactorBase::checkAndNormalizeViewFactor()
188 : {
189 : // collect statistics
190 142 : Real max_sum_deviation = maxDevRowSum();
191 142 : Real max_reciprocity_deviation = maxDevReciprocity();
192 142 : Real max_correction = std::numeric_limits<Real>::lowest();
193 142 : Real min_correction = std::numeric_limits<Real>::max();
194 :
195 142 : if (_print_view_factor_info)
196 0 : _console << "\nSum of all view factors from side i to all other faces before normalization.\n"
197 0 : << std::endl;
198 :
199 : // check view factors
200 142 : if (_print_view_factor_info)
201 0 : for (unsigned int from = 0; from < _n_sides; ++from)
202 0 : _console << "View factors from sideset " << boundaryNames()[from] << " sum to "
203 0 : << viewFactorRowSum(from) << std::endl;
204 :
205 142 : if (max_sum_deviation > _view_factor_tol)
206 0 : mooseError("Maximum deviation of view factor row sum is ",
207 : max_sum_deviation,
208 : " exceeding the set tolerance of ",
209 0 : _view_factor_tol);
210 :
211 : // normalize view factors
212 142 : if (_normalize_view_factor)
213 : {
214 88 : if (_print_view_factor_info)
215 0 : _console << "\nNormalizing view factors.\n" << std::endl;
216 :
217 88 : const auto view_factors_original = _view_factors;
218 :
219 88 : switch (_normalization_method)
220 : {
221 80 : case NormalizationMethod::LAGRANGE_MULTIPLIER:
222 80 : normalizeUsingLagrangeMultiplier();
223 : break;
224 8 : case NormalizationMethod::INVERSE_ROW_SUM:
225 8 : normalizeUsingInverseRowSum();
226 : break;
227 0 : default:
228 0 : mooseError("Invalid normalization method.");
229 : }
230 :
231 : // compute max and min correction for reporting
232 716 : for (const auto i : make_range(_n_sides))
233 5456 : for (const auto j : make_range(_n_sides))
234 : {
235 4828 : const Real correction = _view_factors[i][j] - view_factors_original[i][j];
236 4828 : min_correction = std::min(correction, min_correction);
237 4828 : max_correction = std::max(correction, max_correction);
238 : }
239 88 : }
240 :
241 142 : if (_print_view_factor_info)
242 : {
243 0 : _console << "\nFinal view factors.\n" << std::endl;
244 0 : for (unsigned int from = 0; from < _n_sides; ++from)
245 : {
246 : std::string from_name;
247 0 : for (auto pair : _side_name_index)
248 0 : if (pair.second == from)
249 : from_name = pair.first;
250 0 : auto from_id = _mesh.getBoundaryID(from_name);
251 :
252 0 : for (unsigned int to = 0; to < _n_sides; ++to)
253 : {
254 : std::string to_name;
255 0 : for (auto pair : _side_name_index)
256 0 : if (pair.second == to)
257 : to_name = pair.first;
258 0 : auto to_id = _mesh.getBoundaryID(to_name);
259 0 : _console << from_name << " (" << from_id << ") -> " << to_name << " (" << to_id
260 0 : << ") = " << _view_factors[from][to] << std::endl;
261 : }
262 : }
263 : }
264 :
265 : // check sum of view factors and maximum deviation on reciprocity
266 142 : Real max_sum_deviation_after_normalization = maxDevRowSum();
267 142 : Real max_reciprocity_deviation_after_normalization = maxDevReciprocity();
268 :
269 : // print summary
270 142 : _console << std::endl;
271 142 : _console << COLOR_CYAN << "Summary of the view factor computation"
272 142 : << "\n"
273 142 : << COLOR_DEFAULT << std::endl;
274 142 : if (_normalize_view_factor)
275 88 : _console << "Normalization is performed." << std::endl;
276 : else
277 54 : _console << "Normalization is skipped." << std::endl;
278 142 : _console << std::setw(60) << std::left << "Number of patches: " << _n_sides << std::endl;
279 142 : _console << std::setw(60) << std::left
280 142 : << "Max deviation sum = 1 before normalization: " << max_sum_deviation << std::endl;
281 142 : _console << std::setw(60) << std::left
282 142 : << "Max deviation from reciprocity before normalization: " << max_reciprocity_deviation
283 142 : << std::endl;
284 142 : if (_normalize_view_factor)
285 : {
286 88 : _console << std::setw(60) << std::left << "Maximum correction: " << max_correction << std::endl;
287 88 : _console << std::setw(60) << std::left << "Minimum correction: " << min_correction << std::endl;
288 88 : _console << std::setw(60) << "Max deviation sum = 1 after normalization: "
289 88 : << max_sum_deviation_after_normalization << std::endl;
290 88 : _console << std::setw(60) << std::left << "Max deviation from reciprocity after normalization: "
291 88 : << max_reciprocity_deviation_after_normalization << std::endl;
292 88 : _console << std::endl;
293 : }
294 142 : }
295 :
296 : void
297 80 : ViewFactorBase::normalizeUsingLagrangeMultiplier()
298 : {
299 : // allocate space
300 80 : DenseVector<Real> rhs(_n_sides);
301 80 : DenseVector<Real> lmults(_n_sides);
302 80 : DenseMatrix<Real> matrix(_n_sides, _n_sides);
303 :
304 : // equations for the Lagrange multiplier
305 676 : for (unsigned int i = 0; i < _n_sides; ++i)
306 : {
307 : // set the right hand side
308 596 : rhs(i) = 1 - viewFactorRowSum(i);
309 :
310 : // contribution from the delta_ii element
311 596 : matrix(i, i) = -0.5;
312 :
313 : // contributions from the upper diagonal
314 2648 : for (unsigned int j = i + 1; j < _n_sides; ++j)
315 : {
316 2052 : Real ar = _areas[i] / _areas[j];
317 2052 : Real f = 2 * (1 + ar * ar);
318 2052 : matrix(i, i) += -1 / f;
319 2052 : matrix(i, j) += -1 / f * ar;
320 2052 : rhs(i) += ar * ar / (1 + ar * ar) * devReciprocity(i, j);
321 : }
322 :
323 : // contributions from the lower diagonal
324 2648 : for (unsigned int j = 0; j < i; ++j)
325 : {
326 2052 : Real ar = _areas[j] / _areas[i];
327 2052 : Real f = 2 * (1 + ar * ar);
328 2052 : matrix(i, i) += -1 / f * ar * ar;
329 2052 : matrix(i, j) += -1 / f * ar;
330 2052 : rhs(i) -= ar * devReciprocity(j, i) * (1 - ar * ar / (1 + ar * ar));
331 : }
332 : }
333 :
334 : // solve the linear system
335 80 : matrix.lu_solve(rhs, lmults);
336 :
337 : // perform corrections but store view factors to temporary array
338 : // because we cannot modify _view_factors as it's used in this calc
339 80 : std::vector<std::vector<Real>> vf_temp = _view_factors;
340 676 : for (unsigned int i = 0; i < _n_sides; ++i)
341 5296 : for (unsigned int j = 0; j < _n_sides; ++j)
342 : {
343 : Real correction;
344 4700 : if (i == j)
345 596 : correction = -0.5 * lmults(i);
346 : else
347 : {
348 4104 : Real ar = _areas[i] / _areas[j];
349 4104 : Real f = 2 * (1 + ar * ar);
350 4104 : correction = -(lmults(i) + lmults(j) * ar + 2 * ar * ar * devReciprocity(i, j)) / f;
351 : }
352 :
353 : // update the temporary view factor
354 4700 : vf_temp[i][j] += correction;
355 : }
356 80 : _view_factors = vf_temp;
357 160 : }
358 :
359 : void
360 8 : ViewFactorBase::normalizeUsingInverseRowSum()
361 : {
362 40 : for (unsigned int i = 0; i < _n_sides; ++i)
363 : {
364 32 : const Real row_sum = viewFactorRowSum(i);
365 160 : for (auto & v : _view_factors[i])
366 128 : v /= row_sum;
367 : }
368 8 : }
369 :
370 : unsigned int
371 0 : ViewFactorBase::indexHelper(unsigned int i, unsigned int j) const
372 : {
373 : mooseAssert(i <= j, "indexHelper requires i <= j");
374 0 : if (i == j)
375 0 : return _n_sides + i;
376 0 : unsigned int pos = 2 * _n_sides;
377 0 : for (unsigned int l = 0; l < _n_sides; ++l)
378 0 : for (unsigned int m = l + 1; m < _n_sides; ++m)
379 : {
380 0 : if (l == i && m == j)
381 0 : return pos;
382 : else
383 0 : ++pos;
384 : }
385 0 : mooseError("Should never get here");
386 : return 0;
387 : }
|
