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 "ParameterMeshOptimization.h"
11 :
12 : #include "AddVariableAction.h"
13 : #include "ParameterMesh.h"
14 : #include "OptUtils.h"
15 : #include "libmesh/string_to_enum.h"
16 :
17 : #include "ReadExodusMeshVars.h"
18 :
19 : using namespace libMesh;
20 :
21 : registerMooseObject("OptimizationApp", ParameterMeshOptimization);
22 :
23 : InputParameters
24 251 : ParameterMeshOptimization::validParams()
25 : {
26 251 : InputParameters params = GeneralOptimization::validParams();
27 :
28 251 : params.addClassDescription(
29 : "Computes objective function, gradient and contains reporters for communicating between "
30 : "optimizeSolve and subapps using mesh-based parameter definition.");
31 :
32 502 : params.addRequiredParam<std::vector<FileName>>(
33 : "parameter_meshes", "Exodus file containing meshes describing parameters.");
34 :
35 251 : const auto family = AddVariableAction::getNonlinearVariableFamilies();
36 251 : MultiMooseEnum families(family.getRawNames(), "LAGRANGE");
37 502 : params.addParam<MultiMooseEnum>(
38 : "parameter_families",
39 : families,
40 : "Specifies the family of FE shape functions for each group of parameters. If a single value "
41 : "is "
42 : "specified, then that value is used for all groups of parameters.");
43 251 : const auto order = AddVariableAction::getNonlinearVariableOrders();
44 251 : MultiMooseEnum orders(order.getRawNames(), "FIRST");
45 502 : params.addParam<MultiMooseEnum>(
46 : "parameter_orders",
47 : orders,
48 : "Specifies the order of FE shape functions for each group of parameters. If a single value "
49 : "is "
50 : "specified, then that value is used for all groups of parameters.");
51 :
52 502 : params.addParam<unsigned int>(
53 502 : "num_parameter_times", 1, "The number of time points the parameters represent.");
54 :
55 502 : params.addParam<std::vector<std::string>>(
56 : "initial_condition_mesh_variable",
57 : "Name of variable on parameter mesh to use as initial condition.");
58 502 : params.addParam<std::vector<std::string>>(
59 : "lower_bound_mesh_variable", "Name of variable on parameter mesh to use as lower bound.");
60 502 : params.addParam<std::vector<std::string>>(
61 : "upper_bound_mesh_variable", "Name of variable on parameter mesh to use as upper bound.");
62 502 : params.addParam<std::vector<unsigned int>>(
63 : "exodus_timesteps_for_parameter_mesh_variable",
64 : "Timesteps to read all parameter group bounds and initial conditions from Exodus mesh. The "
65 : "options are to give no timestep, a single timestep or \"num_parameter_times\" timesteps. "
66 : "No timestep results in the final timestep from the mesh being used. A single timestep "
67 : "results in values at that timestep being used for all timesteps. \"num_parameter_times\" "
68 : "timesteps results in values from the mesh at those steps being used. The same timesteps "
69 : "are used for all parameter groups and all meshes, the capability to define different "
70 : "timesteps for different meshes is not supported.");
71 :
72 : // New parameters for multiple regularization types
73 251 : MultiMooseEnum reg_types("L2_GRADIENT");
74 502 : params.addParam<MultiMooseEnum>(
75 : "regularization_types",
76 : reg_types,
77 : "Types of regularization to apply. Multiple types can be specified.");
78 :
79 502 : params.addParam<std::vector<Real>>("regularization_coeffs",
80 : {},
81 : "Coefficients for each regularization type. Must match the "
82 : "number of regularization_types specified.");
83 :
84 502 : params.addParamNamesToGroup("tikhonov_coeff regularization_types regularization_coeffs",
85 : "Regularization");
86 :
87 251 : return params;
88 251 : }
89 :
90 126 : ParameterMeshOptimization::ParameterMeshOptimization(const InputParameters & parameters)
91 : : GeneralOptimization(parameters),
92 252 : _regularization_coeffs(getParam<std::vector<Real>>("regularization_coeffs")),
93 378 : _regularization_types(getParam<MultiMooseEnum>("regularization_types")
94 126 : .getSetValueIDs<ParameterMesh::RegularizationType>())
95 : {
96 : // Validate that regularization coefficients match types
97 126 : if (_regularization_coeffs.size() != _regularization_types.size())
98 2 : paramError("regularization_coeffs",
99 : "Number of regularization coefficients (",
100 : _regularization_coeffs.size(),
101 : ") must match number of regularization types (",
102 : _regularization_types.size(),
103 : ")");
104 124 : }
105 :
106 : std::vector<Real>
107 168 : ParameterMeshOptimization::parseExodusData(const FEType fetype,
108 : const FileName mesh_file_name,
109 : const std::vector<unsigned int> & exodus_timestep,
110 : const std::string & mesh_var_name) const
111 : {
112 : // read data off Exodus mesh
113 334 : ReadExodusMeshVars data_mesh(fetype, mesh_file_name, mesh_var_name);
114 : std::vector<Real> parsed_data;
115 : // read from mesh
116 460 : for (auto const & step : exodus_timestep)
117 : {
118 296 : std::vector<Real> data = data_mesh.getParameterValues(step);
119 294 : parsed_data.insert(parsed_data.end(), data.begin(), data.end());
120 294 : }
121 :
122 164 : return parsed_data;
123 164 : }
124 :
125 : void
126 124 : ParameterMeshOptimization::setICsandBounds()
127 : {
128 496 : if ((isParamValid("num_values_name") || isParamValid("num_values")))
129 0 : paramError("num_values_name or num_values should not be used with ParameterMeshOptimization. "
130 : "Instead the number of dofs is set by the parameter meshes.");
131 :
132 124 : _nvalues.resize(_nparams, 0);
133 : // Fill the mesh information
134 124 : const auto & meshes = getParam<std::vector<FileName>>("parameter_meshes");
135 124 : const auto & families = getParam<MultiMooseEnum>("parameter_families");
136 124 : const auto & orders = getParam<MultiMooseEnum>("parameter_orders");
137 248 : const auto & ntimes = getParam<unsigned int>("num_parameter_times");
138 :
139 : // Fill exodus parameter bounds and IC information
140 : std::vector<std::string> initial_condition_mesh_variable;
141 : std::vector<std::string> lower_bound_mesh_variable;
142 : std::vector<std::string> upper_bound_mesh_variable;
143 248 : if (isParamValid("initial_condition_mesh_variable"))
144 : initial_condition_mesh_variable =
145 156 : getParam<std::vector<std::string>>("initial_condition_mesh_variable");
146 248 : if (isParamValid("lower_bound_mesh_variable"))
147 87 : lower_bound_mesh_variable = getParam<std::vector<std::string>>("lower_bound_mesh_variable");
148 248 : if (isParamValid("upper_bound_mesh_variable"))
149 117 : upper_bound_mesh_variable = getParam<std::vector<std::string>>("upper_bound_mesh_variable");
150 :
151 : std::vector<unsigned int> exodus_timestep;
152 248 : if (isParamValid("exodus_timesteps_for_parameter_mesh_variable"))
153 : exodus_timestep =
154 78 : getParam<std::vector<unsigned int>>("exodus_timesteps_for_parameter_mesh_variable");
155 : else // get last timestep in file
156 98 : exodus_timestep = {std::numeric_limits<unsigned int>::max()};
157 :
158 : // now do a bunch of error checking
159 : // Size checks for data
160 124 : if (meshes.size() != _nparams)
161 1 : paramError("parameter_meshes",
162 : "There must be a mesh associated with each group of parameters.");
163 123 : if (families.size() > 1 && families.size() != _nparams)
164 1 : paramError("parameter_families",
165 : "There must be a family associated with each group of parameters.");
166 122 : if (orders.size() > 1 && orders.size() != _nparams)
167 1 : paramError("parameter_orders",
168 : "There must be an order associated with each group of parameters.");
169 :
170 : // error checking that initial conditions and bounds are only read from a single location
171 346 : if (isParamValid("initial_condition_mesh_variable") && isParamValid("initial_condition"))
172 0 : paramError("initial_condition_mesh_variable",
173 : "Initial conditions for all parameter groups can only be defined by "
174 : "initial_condition_mesh_variable or "
175 : "initial_condition but not both.");
176 300 : else if (isParamValid("lower_bound_mesh_variable") && isParamValid("lower_bounds"))
177 2 : paramError(
178 : "lower_bound_mesh_variable",
179 : "Lower bounds for all parameter groups can only be defined by lower_bound_mesh_variable or "
180 : "lower_bounds but not both.");
181 312 : else if (isParamValid("upper_bound_mesh_variable") && isParamValid("upper_bounds"))
182 0 : paramError(
183 : "upper_bound_mesh_variable",
184 : "Upper bounds for all parameter groups can only be defined by upper_bound_mesh_variable or "
185 : "upper_bounds but not both.");
186 :
187 : // Make sure they did not specify too many timesteps
188 357 : if (isParamValid("exodus_timesteps_for_parameter_mesh_variable") &&
189 197 : (!isParamValid("lower_bound_mesh_variable") + !isParamValid("upper_bound_mesh_variable") +
190 195 : !isParamValid("initial_condition_mesh_variable") ==
191 : 3))
192 2 : paramError("\"exodus_timesteps_for_parameter_mesh_variable\" should only be specified if "
193 : "reading values from a mesh.");
194 117 : else if (exodus_timestep.size() != ntimes && exodus_timestep.size() != 1)
195 2 : paramError("exodus_timesteps_for_parameter_mesh_variable",
196 : "Number of timesteps to read mesh data specified by "
197 : "\"exodus_timesteps_for_parameter_mesh_variable\" incorrect. "
198 : "\"exodus_timesteps_for_parameter_mesh_variable\" can specify a single timestep or "
199 : "\"num_parameter_times\" timesteps.");
200 :
201 115 : _ndof = 0;
202 115 : _parameter_meshes.resize(_nparams);
203 253 : for (const auto & param_id : make_range(_nparams))
204 : {
205 143 : const std::string family = families.size() > 1 ? families[param_id] : families[0];
206 143 : const std::string order = orders.size() > 1 ? orders[param_id] : orders[0];
207 143 : const FEType fetype(Utility::string_to_enum<Order>(order),
208 143 : Utility::string_to_enum<FEFamily>(family));
209 :
210 286 : _parameter_meshes[param_id] = std::make_unique<ParameterMesh>(fetype, meshes[param_id]);
211 143 : _nvalues[param_id] = _parameter_meshes[param_id]->size() * ntimes;
212 143 : _ndof += _nvalues[param_id];
213 :
214 : // read and assign initial conditions
215 : {
216 : std::vector<Real> initial_condition;
217 286 : if (isParamValid("initial_condition_mesh_variable"))
218 219 : initial_condition = parseExodusData(
219 : fetype, meshes[param_id], exodus_timestep, initial_condition_mesh_variable[param_id]);
220 : else
221 208 : initial_condition = parseInputData("initial_condition", 0, param_id);
222 :
223 142 : _parameters[param_id]->assign(initial_condition.begin(), initial_condition.end());
224 142 : }
225 :
226 : // read and assign lower bound
227 : {
228 : std::vector<Real> lower_bound;
229 284 : if (isParamValid("lower_bound_mesh_variable"))
230 109 : lower_bound = parseExodusData(
231 : fetype, meshes[param_id], exodus_timestep, lower_bound_mesh_variable[param_id]);
232 : else
233 309 : lower_bound = parseInputData("lower_bounds", std::numeric_limits<Real>::lowest(), param_id);
234 :
235 138 : _lower_bounds.insert(_lower_bounds.end(), lower_bound.begin(), lower_bound.end());
236 138 : }
237 :
238 : // read and assign upper bound
239 : {
240 : std::vector<Real> upper_bound;
241 276 : if (isParamValid("upper_bound_mesh_variable"))
242 168 : upper_bound = parseExodusData(
243 : fetype, meshes[param_id], exodus_timestep, upper_bound_mesh_variable[param_id]);
244 : else
245 246 : upper_bound = parseInputData("upper_bounds", std::numeric_limits<Real>::max(), param_id);
246 :
247 138 : _upper_bounds.insert(_upper_bounds.end(), upper_bound.begin(), upper_bound.end());
248 138 : }
249 :
250 : // resize gradient vector to be filled later
251 138 : _gradients[param_id]->resize(_nvalues[param_id]);
252 : }
253 110 : }
254 :
255 : Real
256 2272 : ParameterMeshOptimization::computeObjective()
257 : {
258 2272 : Real val = GeneralOptimization::computeObjective();
259 :
260 : // Apply each regularization type with its coefficient
261 4165 : for (const auto reg_idx : index_range(_regularization_types))
262 : {
263 1893 : if (_regularization_coeffs[reg_idx] > 0.0)
264 : {
265 : Real regularization_value = 0.0;
266 :
267 : // Convert MultiMooseEnum to RegularizationType using get() method
268 1893 : ParameterMesh::RegularizationType reg_type = _regularization_types[reg_idx];
269 :
270 3786 : for (const auto & param_id : make_range(_nparams))
271 : {
272 : // Get current parameter values for this group
273 1893 : const auto & param_values = *_parameters[param_id];
274 :
275 : // Compute regularization objective for this type
276 1893 : regularization_value +=
277 1893 : _parameter_meshes[param_id]->computeRegularizationObjective(param_values, reg_type);
278 : }
279 :
280 1893 : val += _regularization_coeffs[reg_idx] * regularization_value;
281 : }
282 : }
283 :
284 2272 : return val;
285 : }
286 :
287 : void
288 1462 : ParameterMeshOptimization::computeGradient(libMesh::PetscVector<Number> & gradient) const
289 : {
290 : // Add regularization gradient contributions to the reporter gradients before base computation
291 2545 : for (const auto reg_idx : index_range(_regularization_types))
292 : {
293 1083 : if (_regularization_coeffs[reg_idx] > 0.0)
294 : {
295 : // Convert MultiMooseEnum to RegularizationType using get() method
296 1083 : ParameterMesh::RegularizationType reg_type = _regularization_types[reg_idx];
297 :
298 2166 : for (const auto & param_id : make_range(_nparams))
299 : {
300 : // Get current parameter values for this group
301 1083 : const auto & param_values = *_parameters[param_id];
302 1083 : auto grad_values = _gradients[param_id];
303 :
304 : // Compute regularization gradient for this type
305 : std::vector<Real> reg_grad =
306 1083 : _parameter_meshes[param_id]->computeRegularizationGradient(param_values, reg_type);
307 :
308 : // Add to gradient with coefficient
309 187320 : for (unsigned int i = 0; i < param_values.size(); ++i)
310 186237 : (*grad_values)[i] += _regularization_coeffs[reg_idx] * reg_grad[i];
311 1083 : }
312 : }
313 : }
314 :
315 : // Now call base class method which includes Tikhonov and copies to PETSc vector
316 1462 : OptimizationReporterBase::computeGradient(gradient);
317 1461 : }
|
