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 338 : ParameterMeshOptimization::validParams()
25 : {
26 338 : InputParameters params = GeneralOptimization::validParams();
27 :
28 338 : params.addClassDescription(
29 : "Computes objective function, gradient and contains reporters for communicating between "
30 : "optimizeSolve and subapps using mesh-based parameter definition.");
31 :
32 676 : params.addRequiredParam<std::vector<FileName>>(
33 : "parameter_meshes", "Exodus file containing meshes describing parameters.");
34 :
35 338 : const auto family = AddVariableAction::getNonlinearVariableFamilies();
36 338 : MultiMooseEnum families(family.getRawNames(), "LAGRANGE");
37 676 : 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 338 : const auto order = AddVariableAction::getNonlinearVariableOrders();
44 338 : MultiMooseEnum orders(order.getRawNames(), "FIRST");
45 676 : 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 676 : params.addParam<unsigned int>(
53 676 : "num_parameter_times", 1, "The number of time points the parameters represent.");
54 :
55 676 : params.addParam<std::vector<std::string>>(
56 : "initial_condition_mesh_variable",
57 : "Name of variable on parameter mesh to use as initial condition.");
58 676 : params.addParam<std::vector<std::string>>(
59 : "lower_bound_mesh_variable", "Name of variable on parameter mesh to use as lower bound.");
60 676 : params.addParam<std::vector<std::string>>(
61 : "upper_bound_mesh_variable", "Name of variable on parameter mesh to use as upper bound.");
62 676 : 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 338 : MultiMooseEnum reg_types("L2_GRADIENT");
74 676 : params.addParam<MultiMooseEnum>(
75 : "regularization_types",
76 : reg_types,
77 : "Types of regularization to apply. Multiple types can be specified.");
78 :
79 676 : 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 676 : params.addParamNamesToGroup("tikhonov_coeff regularization_types regularization_coeffs",
85 : "Regularization");
86 :
87 338 : return params;
88 338 : }
89 :
90 170 : ParameterMeshOptimization::ParameterMeshOptimization(const InputParameters & parameters)
91 : : GeneralOptimization(parameters),
92 340 : _regularization_coeffs(getParam<std::vector<Real>>("regularization_coeffs")),
93 510 : _regularization_types(getParam<MultiMooseEnum>("regularization_types")
94 170 : .getSetValueIDs<ParameterMesh::RegularizationType>())
95 : {
96 : // Validate that regularization coefficients match types
97 170 : if (_regularization_coeffs.size() != _regularization_types.size())
98 4 : 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 166 : }
105 :
106 : std::vector<Real>
107 211 : 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 420 : ReadExodusMeshVars data_mesh(fetype, mesh_file_name, mesh_var_name);
114 : std::vector<Real> parsed_data;
115 : // read from mesh
116 582 : for (auto const & step : exodus_timestep)
117 : {
118 375 : std::vector<Real> data = data_mesh.getParameterValues(step);
119 373 : parsed_data.insert(parsed_data.end(), data.begin(), data.end());
120 373 : }
121 :
122 207 : return parsed_data;
123 207 : }
124 :
125 : void
126 166 : ParameterMeshOptimization::setICsandBounds()
127 : {
128 664 : 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 166 : _nvalues.resize(_nparams, 0);
133 : // Fill the mesh information
134 166 : const auto & meshes = getParam<std::vector<FileName>>("parameter_meshes");
135 166 : const auto & families = getParam<MultiMooseEnum>("parameter_families");
136 166 : const auto & orders = getParam<MultiMooseEnum>("parameter_orders");
137 332 : 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 332 : if (isParamValid("initial_condition_mesh_variable"))
144 : initial_condition_mesh_variable =
145 189 : getParam<std::vector<std::string>>("initial_condition_mesh_variable");
146 332 : if (isParamValid("lower_bound_mesh_variable"))
147 105 : lower_bound_mesh_variable = getParam<std::vector<std::string>>("lower_bound_mesh_variable");
148 332 : if (isParamValid("upper_bound_mesh_variable"))
149 147 : upper_bound_mesh_variable = getParam<std::vector<std::string>>("upper_bound_mesh_variable");
150 :
151 : std::vector<unsigned int> exodus_timestep;
152 332 : if (isParamValid("exodus_timesteps_for_parameter_mesh_variable"))
153 : exodus_timestep =
154 99 : getParam<std::vector<unsigned int>>("exodus_timesteps_for_parameter_mesh_variable");
155 : else // get last timestep in file
156 133 : exodus_timestep = {std::numeric_limits<unsigned int>::max()};
157 :
158 : // now do a bunch of error checking
159 : // Size checks for data
160 166 : if (meshes.size() != _nparams)
161 2 : paramError("parameter_meshes",
162 : "There must be a mesh associated with each group of parameters.");
163 164 : if (families.size() > 1 && families.size() != _nparams)
164 2 : paramError("parameter_families",
165 : "There must be a family associated with each group of parameters.");
166 162 : if (orders.size() > 1 && orders.size() != _nparams)
167 2 : 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 446 : 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 390 : 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 410 : 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 474 : if (isParamValid("exodus_timesteps_for_parameter_mesh_variable") &&
189 257 : (!isParamValid("lower_bound_mesh_variable") + !isParamValid("upper_bound_mesh_variable") +
190 255 : !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 156 : 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 154 : _ndof = 0;
202 154 : _parameter_meshes.resize(_nparams);
203 338 : for (const auto & param_id : make_range(_nparams))
204 : {
205 190 : const std::string family = families.size() > 1 ? families[param_id] : families[0];
206 190 : const std::string order = orders.size() > 1 ? orders[param_id] : orders[0];
207 190 : const FEType fetype(Utility::string_to_enum<Order>(order),
208 190 : Utility::string_to_enum<FEFamily>(family));
209 :
210 380 : _parameter_meshes[param_id] = std::make_unique<ParameterMesh>(fetype, meshes[param_id]);
211 190 : _nvalues[param_id] = _parameter_meshes[param_id]->size() * ntimes;
212 190 : _ndof += _nvalues[param_id];
213 :
214 : // read and assign initial conditions
215 : {
216 : std::vector<Real> initial_condition;
217 380 : if (isParamValid("initial_condition_mesh_variable"))
218 270 : initial_condition = parseExodusData(
219 : fetype, meshes[param_id], exodus_timestep, initial_condition_mesh_variable[param_id]);
220 : else
221 296 : initial_condition = parseInputData("initial_condition", 0, param_id);
222 :
223 188 : _parameters[param_id]->assign(initial_condition.begin(), initial_condition.end());
224 188 : }
225 :
226 : // read and assign lower bound
227 : {
228 : std::vector<Real> lower_bound;
229 376 : if (isParamValid("lower_bound_mesh_variable"))
230 139 : lower_bound = parseExodusData(
231 : fetype, meshes[param_id], exodus_timestep, lower_bound_mesh_variable[param_id]);
232 : else
233 417 : lower_bound = parseInputData("lower_bounds", std::numeric_limits<Real>::lowest(), param_id);
234 :
235 184 : _lower_bounds.insert(_lower_bounds.end(), lower_bound.begin(), lower_bound.end());
236 184 : }
237 :
238 : // read and assign upper bound
239 : {
240 : std::vector<Real> upper_bound;
241 368 : if (isParamValid("upper_bound_mesh_variable"))
242 216 : upper_bound = parseExodusData(
243 : fetype, meshes[param_id], exodus_timestep, upper_bound_mesh_variable[param_id]);
244 : else
245 336 : upper_bound = parseInputData("upper_bounds", std::numeric_limits<Real>::max(), param_id);
246 :
247 184 : _upper_bounds.insert(_upper_bounds.end(), upper_bound.begin(), upper_bound.end());
248 184 : }
249 :
250 : // resize gradient vector to be filled later
251 184 : _gradients[param_id]->resize(_nvalues[param_id]);
252 : }
253 148 : }
254 :
255 : Real
256 2171 : ParameterMeshOptimization::computeObjective()
257 : {
258 2171 : Real val = GeneralOptimization::computeObjective();
259 :
260 : // Apply each regularization type with its coefficient
261 3773 : for (const auto reg_idx : index_range(_regularization_types))
262 : {
263 1602 : if (_regularization_coeffs[reg_idx] > 0.0)
264 : {
265 : Real regularization_value = 0.0;
266 :
267 : // Convert MultiMooseEnum to RegularizationType using get() method
268 1602 : ParameterMesh::RegularizationType reg_type = _regularization_types[reg_idx];
269 :
270 3204 : for (const auto & param_id : make_range(_nparams))
271 : {
272 : // Get current parameter values for this group
273 1602 : const auto & param_values = *_parameters[param_id];
274 :
275 : // Compute regularization objective for this type
276 1602 : regularization_value +=
277 1602 : _parameter_meshes[param_id]->computeRegularizationObjective(param_values, reg_type);
278 : }
279 :
280 1602 : val += _regularization_coeffs[reg_idx] * regularization_value;
281 : }
282 : }
283 :
284 2171 : return val;
285 : }
286 :
287 : void
288 2171 : ParameterMeshOptimization::computeGradient(libMesh::PetscVector<Number> & gradient) const
289 : {
290 : // Add regularization gradient contributions to the reporter gradients before base computation
291 3773 : for (const auto reg_idx : index_range(_regularization_types))
292 : {
293 1602 : if (_regularization_coeffs[reg_idx] > 0.0)
294 : {
295 : // Convert MultiMooseEnum to RegularizationType using get() method
296 1602 : ParameterMesh::RegularizationType reg_type = _regularization_types[reg_idx];
297 :
298 3204 : for (const auto & param_id : make_range(_nparams))
299 : {
300 : // Get current parameter values for this group
301 1602 : const auto & param_values = *_parameters[param_id];
302 1602 : auto grad_values = _gradients[param_id];
303 :
304 : // Compute regularization gradient for this type
305 : std::vector<Real> reg_grad =
306 1602 : _parameter_meshes[param_id]->computeRegularizationGradient(param_values, reg_type);
307 :
308 : // Add to gradient with coefficient
309 280350 : for (unsigned int i = 0; i < param_values.size(); ++i)
310 278748 : (*grad_values)[i] += _regularization_coeffs[reg_idx] * reg_grad[i];
311 1602 : }
312 : }
313 : }
314 :
315 : // Now call base class method which includes Tikhonov and copies to PETSc vector
316 2171 : OptimizationReporterBase::computeGradient(gradient);
317 2169 : }
|
