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 "PODReducedBasisSurrogate.h"
11 :
12 : registerMooseObject("StochasticToolsApp", PODReducedBasisSurrogate);
13 :
14 : InputParameters
15 132 : PODReducedBasisSurrogate::validParams()
16 : {
17 132 : InputParameters params = SurrogateModel::validParams();
18 132 : params.addClassDescription("Evaluates POD-RB surrogate model with reduced operators "
19 : "computed from PODReducedBasisTrainer.");
20 264 : params.addParam<std::vector<std::string>>("change_rank",
21 264 : std::vector<std::string>(0),
22 : "Names of variables whose rank should be changed.");
23 396 : params.addParam<std::vector<unsigned int>>(
24 : "new_ranks",
25 264 : std::vector<unsigned int>(0),
26 : "The new ranks that each variable in 'change_rank' shall have.");
27 264 : params.addParam<Real>("penalty", 1e5, "The penalty parameter for Dirichlet BCs.");
28 132 : return params;
29 0 : }
30 :
31 68 : PODReducedBasisSurrogate::PODReducedBasisSurrogate(const InputParameters & parameters)
32 : : SurrogateModel(parameters),
33 68 : _change_rank(getParam<std::vector<std::string>>("change_rank")),
34 204 : _new_ranks(getParam<std::vector<unsigned int>>("new_ranks")),
35 136 : _var_names(getModelData<std::vector<std::string>>("_var_names")),
36 136 : _tag_types(getModelData<std::vector<std::string>>("_tag_types")),
37 136 : _base(getModelData<std::vector<std::vector<DenseVector<Real>>>>("_base")),
38 136 : _red_operators(getModelData<std::vector<DenseMatrix<Real>>>("_red_operators")),
39 136 : _penalty(getParam<Real>("penalty")),
40 204 : _initialized(false)
41 : {
42 68 : if (_change_rank.size() != _new_ranks.size())
43 4 : paramError("new_ranks",
44 : "The size of 'new_ranks' is not equal to the ",
45 : "size of 'change_rank' ",
46 : _new_ranks.size(),
47 : " != ",
48 : _change_rank.size());
49 :
50 80 : for (unsigned int var_i = 0; var_i < _new_ranks.size(); ++var_i)
51 16 : if (_new_ranks[var_i] == 0)
52 0 : paramError("new_ranks", "The values should be greater than 0!");
53 64 : }
54 :
55 : void
56 320 : PODReducedBasisSurrogate::evaluateSolution(const std::vector<Real> & params)
57 : {
58 320 : if (!_initialized)
59 : {
60 : // The containers are initialized (if needed).
61 64 : initializeReducedSystem();
62 64 : initializeApproximateSolution();
63 64 : _initialized = true;
64 : }
65 :
66 : // Assembling and solving the reduced equation system.
67 320 : solveReducedSystem(params);
68 :
69 : // Reconstructing the approximate solutions for every variable.
70 320 : reconstructApproximateSolution();
71 320 : }
72 :
73 : void
74 0 : PODReducedBasisSurrogate::evaluateSolution(const std::vector<Real> & params,
75 : DenseVector<Real> & inp_vector,
76 : std::string var_name)
77 : {
78 0 : if (!_initialized)
79 : {
80 : // The containers are initialized (if needed).
81 0 : initializeReducedSystem();
82 0 : _initialized = true;
83 : }
84 :
85 : // Assembling and solving the reduced equation system.
86 0 : solveReducedSystem(params);
87 :
88 : // Reconstructing the approximate solutions for every variable.
89 0 : reconstructApproximateSolution(inp_vector, var_name);
90 0 : }
91 :
92 : void
93 64 : PODReducedBasisSurrogate::initializeReducedSystem()
94 : {
95 : // Storing important indices for the assembly loops.
96 64 : _final_ranks.resize(_var_names.size());
97 64 : _comulative_ranks.resize(_var_names.size());
98 : unsigned int sum_ranks = 0;
99 :
100 : // Checking if the user wants to overwrite the original ranks for the
101 : // variables.
102 128 : for (unsigned int var_i = 0; var_i < _var_names.size(); ++var_i)
103 : {
104 64 : _final_ranks[var_i] = _base[var_i].size();
105 80 : for (unsigned int var_j = 0; var_j < _change_rank.size(); ++var_j)
106 : {
107 16 : if (_change_rank[var_j] == _var_names[var_i])
108 : {
109 16 : if (_new_ranks[var_j] > _base[var_i].size())
110 : {
111 0 : mooseWarning("The specified new rank (",
112 : _new_ranks[var_j],
113 : ") for variable '",
114 : _var_names[var_i],
115 : "' is higher than the original rank (",
116 0 : _base[var_i].size(),
117 : ")! Switched to original rank.");
118 0 : break;
119 : }
120 :
121 16 : _final_ranks[var_i] = _new_ranks[var_j];
122 : }
123 : }
124 64 : sum_ranks += _final_ranks[var_i];
125 64 : _comulative_ranks[var_i] = sum_ranks;
126 : }
127 :
128 : // Resizing containers to match the newly prescribed ranks.
129 64 : _sys_mx = DenseMatrix<Real>(sum_ranks, sum_ranks);
130 64 : _rhs = DenseVector<Real>(sum_ranks);
131 64 : _coeffs = DenseVector<Real>(sum_ranks);
132 64 : }
133 :
134 : void
135 64 : PODReducedBasisSurrogate::initializeApproximateSolution()
136 : {
137 64 : _approx_solution.resize(_var_names.size());
138 128 : for (unsigned int var_i = 0; var_i < _var_names.size(); var_i++)
139 128 : _approx_solution[var_i] = DenseVector<Real>(_base[var_i][0].size());
140 64 : }
141 :
142 : void
143 320 : PODReducedBasisSurrogate::solveReducedSystem(const std::vector<Real> & params)
144 : {
145 : // Cleaning the containers of the system matrix and right hand side.
146 320 : _sys_mx.zero();
147 : _rhs.zero();
148 :
149 : // The assumption here is that the reduced operators in the trainer were
150 : // assembled in the order of the parameters. Also, if the number of
151 : // parameters is fewer than the number of operators, the operator will
152 : // just be added without scaling.
153 1640 : for (unsigned int i = 0; i < _red_operators.size(); ++i)
154 : {
155 : unsigned int row_start = 0;
156 :
157 : // Checking if the reduced operator corresponds to a Dirichlet BC, if
158 : // yes introduce the penalty factor.
159 : Real factor = 1.0;
160 1320 : if (_tag_types[i] == "op_dir" || _tag_types[i] == "src_dir")
161 200 : factor = _penalty;
162 :
163 : // If the user decreased the rank of the reduced bases manually, some parts
164 : // of the initial reduced operators have to be omited.
165 2640 : for (unsigned int var_i = 0; var_i < _var_names.size(); ++var_i)
166 : {
167 5840 : for (unsigned int row_i = row_start; row_i < _comulative_ranks[var_i]; row_i++)
168 : {
169 4520 : if (_tag_types[i] == "op" || _tag_types[i] == "op_dir")
170 : {
171 :
172 : unsigned int col_start = 0;
173 :
174 5720 : for (unsigned int var_j = 0; var_j < _var_names.size(); ++var_j)
175 : {
176 14120 : for (unsigned int col_i = col_start; col_i < _comulative_ranks[var_j]; col_i++)
177 : {
178 11260 : if (i < params.size())
179 8760 : _sys_mx(row_i, col_i) += params[i] * factor * _red_operators[i](row_i, col_i);
180 : else
181 2500 : _sys_mx(row_i, col_i) += factor * _red_operators[i](row_i, col_i);
182 : }
183 :
184 : col_start = _comulative_ranks[var_j];
185 : }
186 : }
187 : else
188 : {
189 1660 : if (i < params.size())
190 1660 : _rhs(row_i) -= params[i] * factor * _red_operators[i](row_i, 0);
191 : else
192 0 : _rhs(row_i) -= factor * _red_operators[i](row_i, 0);
193 : }
194 4520 : row_start = _comulative_ranks[var_i];
195 : }
196 : }
197 : }
198 :
199 : // Solving the reduced system.
200 320 : _sys_mx.lu_solve(_rhs, _coeffs);
201 320 : }
202 :
203 : void
204 320 : PODReducedBasisSurrogate::reconstructApproximateSolution()
205 : {
206 : unsigned int counter = 0;
207 640 : for (unsigned int var_i = 0; var_i < _var_names.size(); var_i++)
208 : {
209 : _approx_solution[var_i].zero();
210 :
211 : // This also takes into account the potential truncation of the bases by
212 : // the user.
213 1440 : for (unsigned int base_i = 0; base_i < _final_ranks[var_i]; ++base_i)
214 : {
215 25200 : for (unsigned int dof_i = 0; dof_i < _base[var_i][base_i].size(); ++dof_i)
216 24080 : _approx_solution[var_i](dof_i) += _coeffs(counter) * _base[var_i][base_i](dof_i);
217 1120 : counter++;
218 : }
219 : }
220 320 : }
221 :
222 : void
223 0 : PODReducedBasisSurrogate::reconstructApproximateSolution(DenseVector<Real> & inp_vector,
224 : std::string var_name)
225 : {
226 0 : auto it = std::find(_var_names.begin(), _var_names.end(), var_name);
227 0 : if (it == _var_names.end())
228 0 : mooseError("Variable '", var_name, "' does not exist in the POD-RB surrogate!");
229 :
230 0 : unsigned int var_i = std::distance(_var_names.begin(), it);
231 :
232 0 : if (inp_vector.size() != _base[var_i][0].size())
233 0 : mooseError("The size of the input vector (",
234 0 : inp_vector.size(),
235 : ") for variable '",
236 : var_name,
237 : "' does not match the size of the stored base vector (",
238 0 : _base[var_i][0].size(),
239 : ") in POD-RB surrogate!");
240 :
241 : inp_vector.zero();
242 :
243 : unsigned int base_begin = 0;
244 0 : if (var_i != 0)
245 0 : base_begin = _comulative_ranks[var_i - 1];
246 :
247 0 : for (unsigned int base_i = 0; base_i < _final_ranks[var_i]; ++base_i)
248 : {
249 0 : for (unsigned int dof_i = 0; dof_i < inp_vector.size(); ++dof_i)
250 0 : inp_vector(dof_i) += _coeffs(base_i + base_begin) * _base[var_i][base_i](dof_i);
251 : }
252 0 : }
253 :
254 : Real
255 320 : PODReducedBasisSurrogate::getNodalQoI(std::string var_name, unsigned int qoi_type) const
256 : {
257 : Real val = 0.0;
258 :
259 320 : auto it = std::find(_var_names.begin(), _var_names.end(), var_name);
260 320 : if (it == _var_names.end())
261 0 : mooseError("Variable '", var_name, "' not found!");
262 :
263 320 : switch (qoi_type)
264 : {
265 300 : case 0:
266 300 : val = _approx_solution[it - _var_names.begin()].max();
267 300 : break;
268 :
269 0 : case 1:
270 0 : val = _approx_solution[it - _var_names.begin()].min();
271 0 : break;
272 :
273 0 : case 2:
274 0 : val = _approx_solution[it - _var_names.begin()].l1_norm();
275 : break;
276 :
277 20 : case 3:
278 20 : val = _approx_solution[it - _var_names.begin()].l2_norm();
279 : break;
280 :
281 0 : case 4:
282 0 : val = _approx_solution[it - _var_names.begin()].linfty_norm();
283 : break;
284 : }
285 :
286 320 : return (val);
287 : }
|
