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 "PODMapping.h"
11 : #include <slepcsvd.h>
12 : #include "MooseTypes.h"
13 : #include <petscdmda.h>
14 :
15 : registerMooseObject("StochasticToolsApp", PODMapping);
16 :
17 : InputParameters
18 340 : PODMapping::validParams()
19 : {
20 340 : InputParameters params = VariableMappingBase::validParams();
21 340 : params.addClassDescription("Class which provides a Proper Orthogonal Decomposition-based mapping "
22 : "between full-order and reduced-order spaces.");
23 680 : params.addParam<UserObjectName>(
24 : "solution_storage", "The name of the storage reporter where the snapshots are located.");
25 680 : params.addParam<std::vector<dof_id_type>>(
26 : "num_modes_to_compute",
27 : "The number of modes that this object should compute. "
28 : "Modes with 0 eigenvalues are filtered out, so the real number of modes "
29 : "might be lower than this. This is also used for setting the "
30 : "subspace sizes for distributed singular value solves. By default, the subspace used for the "
31 : "SVD is twice as big as the number of requested vectors. For more information see the SLEPc "
32 : "manual. If not specified, only one mode is computed per variable.");
33 680 : params.addParam<std::vector<Real>>(
34 : "energy_threshold",
35 340 : std::vector<Real>(),
36 : "The energy threshold for the automatic truncation of the number of modes. In general, the "
37 : "lower this number is the more information is retained about the system by keeping more POD "
38 : "modes.");
39 680 : params.addParam<std::string>("extra_slepc_options",
40 : "",
41 : "Additional options for the singular/eigenvalue solvers in SLEPc.");
42 340 : return params;
43 0 : }
44 :
45 170 : PODMapping::PODMapping(const InputParameters & parameters)
46 : : VariableMappingBase(parameters),
47 : UserObjectInterface(this),
48 510 : _num_modes(isParamValid("num_modes_to_compute")
49 170 : ? getParam<std::vector<dof_id_type>>("num_modes_to_compute")
50 170 : : std::vector<dof_id_type>(_variable_names.size(), 1)),
51 340 : _energy_threshold(getParam<std::vector<Real>>("energy_threshold")),
52 340 : _left_basis_functions(declareModelData<std::map<VariableName, std::vector<DenseVector<Real>>>>(
53 : "left_basis_functions")),
54 340 : _right_basis_functions(declareModelData<std::map<VariableName, std::vector<DenseVector<Real>>>>(
55 : "right_basis_functions")),
56 170 : _singular_values(
57 170 : declareModelData<std::map<VariableName, std::vector<Real>>>("singular_values")),
58 340 : _extra_slepc_options(getParam<std::string>("extra_slepc_options")),
59 170 : _parallel_storage(isParamValid("solution_storage")
60 260 : ? &getUserObject<ParallelSolutionStorage>("solution_storage")
61 : : nullptr),
62 340 : _pod(StochasticTools::POD(_parallel_storage, _extra_slepc_options, _communicator))
63 : {
64 340 : if (!isParamValid("filename"))
65 : {
66 90 : if (_num_modes.size() != _variable_names.size())
67 0 : paramError("num_modes", "The number of modes should be defined for each variable!");
68 :
69 220 : for (const auto & mode : _num_modes)
70 130 : if (!mode)
71 0 : paramError("num_modes", "The number of modes should always be a positive integer!");
72 :
73 90 : if (_energy_threshold.size())
74 : {
75 0 : if (_energy_threshold.size() != _variable_names.size())
76 0 : paramError("energy_threshold",
77 : "The energy thresholds should be defined for each variable!");
78 :
79 0 : for (const auto & threshold : _energy_threshold)
80 0 : if (threshold < 0 || threshold >= 1)
81 0 : paramError("energy_threshold",
82 : "The energy thresholds should always be in the [0,1) range!");
83 : }
84 :
85 : #if PETSC_VERSION_LESS_THAN(3, 14, 0)
86 : mooseError("PODMapping is not supported with PETSc version below 3.14!");
87 : #else
88 220 : for (const auto & vname : _variable_names)
89 : {
90 130 : _singular_values.emplace(vname, std::vector<Real>());
91 130 : _left_basis_functions.emplace(vname, std::vector<DenseVector<Real>>());
92 130 : _right_basis_functions.emplace(vname, std::vector<DenseVector<Real>>());
93 130 : _mapping_ready_to_use.emplace(vname, false);
94 : }
95 : #endif
96 : }
97 170 : }
98 :
99 : void
100 130 : PODMapping::buildMapping(const VariableName &
101 : #if !PETSC_VERSION_LESS_THAN(3, 14, 0)
102 : vname
103 : #endif
104 : )
105 : {
106 : #if !PETSC_VERSION_LESS_THAN(3, 14, 0)
107 130 : if (!_parallel_storage)
108 0 : paramError(
109 : "solution_storage",
110 : "The parallel storage reporter is not supplied! We cannot build a mapping without data!");
111 :
112 130 : if (_mapping_ready_to_use.find(vname) != _mapping_ready_to_use.end() &&
113 130 : !_mapping_ready_to_use[vname])
114 : {
115 130 : auto it = std::find(_variable_names.begin(), _variable_names.end(), vname);
116 : mooseAssert(it != _variable_names.end(), "Variable " + vname + " is not in PODMapping!");
117 : unsigned int var_i = std::distance(_variable_names.begin(), it);
118 : // Clear storage containers for the basis functions and singular values.
119 130 : _left_basis_functions[vname].clear();
120 130 : _right_basis_functions[vname].clear();
121 130 : _singular_values[vname].clear();
122 : // Find the number of modes that we would want to compute
123 130 : std::size_t num_modes_compute = (std::size_t)_num_modes[var_i];
124 : // Find the first element that satisfies the threshold
125 130 : const Real threshold = (_energy_threshold.empty() ? 0.0 : _energy_threshold[var_i]) +
126 130 : std::numeric_limits<Real>::epsilon();
127 : // Use POD class to compute for a variable
128 130 : _pod.computePOD(vname,
129 130 : _left_basis_functions[vname],
130 130 : _right_basis_functions[vname],
131 130 : _singular_values[vname],
132 : num_modes_compute,
133 : threshold);
134 130 : _mapping_ready_to_use[vname] = true;
135 : }
136 :
137 : #endif
138 130 : }
139 :
140 : void
141 320 : PODMapping::map(const VariableName & vname,
142 : const unsigned int global_sample_i,
143 : std::vector<Real> & reduced_order_vector) const
144 : {
145 : mooseAssert(_parallel_storage, "We need the parallel solution storage for this operation.");
146 : mooseAssert(_left_basis_functions.find(vname) != _left_basis_functions.end(),
147 : "The bases for the requested variable are not available!");
148 320 : checkIfReadyToUse(vname);
149 :
150 : // This takes the 0th snapshot because we only support steady-state simulations
151 : // at the moment.
152 320 : const auto & snapshot = _parallel_storage->getGlobalSample(global_sample_i, vname)[0];
153 :
154 320 : const auto & bases = _left_basis_functions[vname];
155 :
156 : reduced_order_vector.clear();
157 320 : reduced_order_vector.resize(bases.size());
158 880 : for (auto base_i : index_range(bases))
159 560 : reduced_order_vector[base_i] = bases[base_i].dot(snapshot);
160 320 : }
161 :
162 : void
163 160 : PODMapping::map(const VariableName & vname,
164 : const DenseVector<Real> & full_order_vector,
165 : std::vector<Real> & reduced_order_vector) const
166 : {
167 160 : checkIfReadyToUse(vname);
168 :
169 : mooseAssert(_left_basis_functions.find(vname) != _left_basis_functions.end(),
170 : "The bases for the requested variable are not available!");
171 :
172 160 : const auto & bases = _left_basis_functions[vname];
173 :
174 : reduced_order_vector.clear();
175 160 : reduced_order_vector.resize(bases.size());
176 400 : for (auto base_i : index_range(bases))
177 240 : reduced_order_vector[base_i] = bases[base_i].dot(full_order_vector);
178 160 : }
179 :
180 : void
181 120 : PODMapping::inverse_map(const VariableName & vname,
182 : const std::vector<Real> & reduced_order_vector,
183 : DenseVector<Real> & full_order_vector) const
184 : {
185 : mooseAssert(std::find(_variable_names.begin(), _variable_names.end(), vname) !=
186 : _variable_names.end(),
187 : "Variable " + vname + " is not in PODMapping!");
188 :
189 120 : checkIfReadyToUse(vname);
190 :
191 120 : if (reduced_order_vector.size() != _left_basis_functions[vname].size())
192 0 : mooseError("The number of supplied reduced-order coefficients (",
193 0 : reduced_order_vector.size(),
194 : ") is not the same as the number of basis functions (",
195 0 : _left_basis_functions[vname].size(),
196 : ") for variable ",
197 : vname,
198 : "!");
199 : // This zeros the DenseVector too
200 120 : full_order_vector.resize(_left_basis_functions[vname][0].size());
201 :
202 360 : for (auto base_i : index_range(reduced_order_vector))
203 2880 : for (unsigned int dof_i = 0; dof_i < _left_basis_functions[vname][base_i].size(); ++dof_i)
204 2640 : full_order_vector(dof_i) +=
205 2640 : reduced_order_vector[base_i] * _left_basis_functions[vname][base_i](dof_i);
206 120 : }
207 :
208 : const DenseVector<Real> &
209 0 : PODMapping::leftBasisFunction(const VariableName & vname, const unsigned int base_i)
210 : {
211 : mooseAssert(std::find(_variable_names.begin(), _variable_names.end(), vname) !=
212 : _variable_names.end(),
213 : "Variable " + vname + " is not in PODMapping!");
214 :
215 0 : checkIfReadyToUse(vname);
216 :
217 : mooseAssert(base_i < _left_basis_functions[vname].size(),
218 : "The POD for " + vname + " only has " +
219 : std::to_string(_left_basis_functions[vname].size()) + " left modes!");
220 :
221 0 : return _left_basis_functions[vname][base_i];
222 : }
223 :
224 : const DenseVector<Real> &
225 0 : PODMapping::rightBasisFunction(const VariableName & vname, const unsigned int base_i)
226 : {
227 : mooseAssert(std::find(_variable_names.begin(), _variable_names.end(), vname) !=
228 : _variable_names.end(),
229 : "Variable " + vname + " is not in PODMapping!");
230 :
231 0 : checkIfReadyToUse(vname);
232 :
233 : mooseAssert(base_i < _right_basis_functions[vname].size(),
234 : "The POD for " + vname + " only has " +
235 : std::to_string(_right_basis_functions[vname].size()) + " right modes!");
236 :
237 0 : return _right_basis_functions[vname][base_i];
238 : }
239 :
240 : const std::vector<DenseVector<Real>> &
241 50 : PODMapping::leftBasis(const VariableName & vname)
242 : {
243 50 : checkIfReadyToUse(vname);
244 50 : if (_left_basis_functions.find(vname) == _left_basis_functions.end())
245 0 : mooseError("We are trying to access container for variable '",
246 : vname,
247 : "' but we don't have it in the POD mapping!");
248 50 : return _left_basis_functions[vname];
249 : }
250 :
251 : const std::vector<DenseVector<Real>> &
252 50 : PODMapping::rightBasis(const VariableName & vname)
253 : {
254 50 : checkIfReadyToUse(vname);
255 50 : if (_right_basis_functions.find(vname) == _right_basis_functions.end())
256 0 : mooseError("We are trying to access container for variable '",
257 : vname,
258 : "' but we don't have it in the POD mapping!");
259 50 : return _right_basis_functions[vname];
260 : }
261 :
262 : const std::vector<Real> &
263 50 : PODMapping::singularValues(const VariableName & vname)
264 : {
265 50 : checkIfReadyToUse(vname);
266 50 : if (_singular_values.find(vname) == _singular_values.end())
267 0 : mooseError("We are trying to access container for variable '",
268 : vname,
269 : "' but we don't have it in the POD mapping!");
270 50 : return _singular_values[vname];
271 : }
|
