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 168 : PODMapping::validParams()
19 : {
20 168 : InputParameters params = VariableMappingBase::validParams();
21 168 : params.addClassDescription("Class which provides a Proper Orthogonal Decomposition-based mapping "
22 : "between full-order and reduced-order spaces.");
23 336 : params.addParam<UserObjectName>(
24 : "solution_storage", "The name of the storage reporter where the snapshots are located.");
25 336 : 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 168 : params.addParam<std::vector<Real>>(
34 : "energy_threshold",
35 168 : 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 336 : params.addParam<std::string>("extra_slepc_options",
40 : "",
41 : "Additional options for the singular/eigenvalue solvers in SLEPc.");
42 168 : return params;
43 0 : }
44 :
45 84 : PODMapping::PODMapping(const InputParameters & parameters)
46 : : VariableMappingBase(parameters),
47 : UserObjectInterface(this),
48 252 : _num_modes(isParamValid("num_modes_to_compute")
49 84 : ? getParam<std::vector<dof_id_type>>("num_modes_to_compute")
50 84 : : std::vector<dof_id_type>(_variable_names.size(), 1)),
51 168 : _energy_threshold(getParam<std::vector<Real>>("energy_threshold")),
52 168 : _left_basis_functions(declareModelData<std::map<VariableName, std::vector<DenseVector<Real>>>>(
53 : "left_basis_functions")),
54 168 : _right_basis_functions(declareModelData<std::map<VariableName, std::vector<DenseVector<Real>>>>(
55 : "right_basis_functions")),
56 84 : _singular_values(
57 84 : declareModelData<std::map<VariableName, std::vector<Real>>>("singular_values")),
58 168 : _extra_slepc_options(getParam<std::string>("extra_slepc_options")),
59 84 : _parallel_storage(isParamValid("solution_storage")
60 128 : ? &getUserObject<ParallelSolutionStorage>("solution_storage")
61 : : nullptr),
62 168 : _pod(StochasticTools::POD(_parallel_storage, _extra_slepc_options, _communicator))
63 : {
64 168 : if (!isParamValid("filename"))
65 : {
66 44 : if (_num_modes.size() != _variable_names.size())
67 0 : paramError("num_modes", "The number of modes should be defined for each variable!");
68 :
69 108 : for (const auto & mode : _num_modes)
70 64 : if (!mode)
71 0 : paramError("num_modes", "The number of modes should always be a positive integer!");
72 :
73 44 : 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 108 : for (const auto & vname : _variable_names)
89 : {
90 64 : _singular_values.emplace(vname, std::vector<Real>());
91 64 : _left_basis_functions.emplace(vname, std::vector<DenseVector<Real>>());
92 64 : _right_basis_functions.emplace(vname, std::vector<DenseVector<Real>>());
93 64 : _mapping_ready_to_use.emplace(vname, false);
94 : }
95 : #endif
96 : }
97 84 : }
98 :
99 : void
100 64 : 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 64 : 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 64 : if (_mapping_ready_to_use.find(vname) != _mapping_ready_to_use.end() &&
113 64 : !_mapping_ready_to_use[vname])
114 : {
115 64 : 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 64 : _left_basis_functions[vname].clear();
120 64 : _right_basis_functions[vname].clear();
121 64 : _singular_values[vname].clear();
122 : // Find the number of modes that we would want to compute
123 64 : std::size_t num_modes_compute = (std::size_t)_num_modes[var_i];
124 : // Find the first element that satisfies the threshold
125 64 : const Real threshold = (_energy_threshold.empty() ? 0.0 : _energy_threshold[var_i]) +
126 64 : std::numeric_limits<Real>::epsilon();
127 : // Use POD class to compute for a variable
128 64 : _pod.computePOD(vname,
129 64 : _left_basis_functions[vname],
130 64 : _right_basis_functions[vname],
131 64 : _singular_values[vname],
132 : num_modes_compute,
133 : threshold);
134 64 : _mapping_ready_to_use[vname] = true;
135 : }
136 :
137 : #endif
138 64 : }
139 :
140 : void
141 160 : 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 160 : checkIfReadyToUse(vname);
149 :
150 : // This takes the 0th snapshot because we only support steady-state simulations
151 : // at the moment.
152 160 : const auto & snapshot = _parallel_storage->getGlobalSample(global_sample_i, vname)[0];
153 :
154 160 : const auto & bases = _left_basis_functions[vname];
155 :
156 160 : reduced_order_vector.clear();
157 160 : reduced_order_vector.resize(bases.size());
158 440 : for (auto base_i : index_range(bases))
159 280 : reduced_order_vector[base_i] = bases[base_i].dot(snapshot);
160 160 : }
161 :
162 : void
163 80 : PODMapping::map(const VariableName & vname,
164 : const DenseVector<Real> & full_order_vector,
165 : std::vector<Real> & reduced_order_vector) const
166 : {
167 80 : 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 80 : const auto & bases = _left_basis_functions[vname];
173 :
174 80 : reduced_order_vector.clear();
175 80 : reduced_order_vector.resize(bases.size());
176 200 : for (auto base_i : index_range(bases))
177 120 : reduced_order_vector[base_i] = bases[base_i].dot(full_order_vector);
178 80 : }
179 :
180 : void
181 60 : 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 60 : checkIfReadyToUse(vname);
190 :
191 60 : 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 60 : full_order_vector.resize(_left_basis_functions[vname][0].size());
201 :
202 180 : for (auto base_i : index_range(reduced_order_vector))
203 1440 : for (unsigned int dof_i = 0; dof_i < _left_basis_functions[vname][base_i].size(); ++dof_i)
204 1320 : full_order_vector(dof_i) +=
205 1320 : reduced_order_vector[base_i] * _left_basis_functions[vname][base_i](dof_i);
206 60 : }
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 24 : PODMapping::leftBasis(const VariableName & vname)
242 : {
243 24 : checkIfReadyToUse(vname);
244 24 : 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 24 : return _left_basis_functions[vname];
249 : }
250 :
251 : const std::vector<DenseVector<Real>> &
252 24 : PODMapping::rightBasis(const VariableName & vname)
253 : {
254 24 : checkIfReadyToUse(vname);
255 24 : 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 24 : return _right_basis_functions[vname];
260 : }
261 :
262 : const std::vector<Real> &
263 24 : PODMapping::singularValues(const VariableName & vname)
264 : {
265 24 : checkIfReadyToUse(vname);
266 24 : 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 24 : return _singular_values[vname];
271 : }
|
