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 "DensityUpdateTwoConstraints.h"
11 : #include "Transient.h"
12 : #include <algorithm>
13 :
14 : registerMooseObject("OptimizationApp", DensityUpdateTwoConstraints);
15 :
16 : InputParameters
17 38 : DensityUpdateTwoConstraints::validParams()
18 : {
19 38 : InputParameters params = ElementUserObject::validParams();
20 38 : params.addClassDescription(
21 : "Compute updated densities based on sensitivities using an optimality criteria method to "
22 : "keep the volume and cost constraints satisified.");
23 76 : params.addParam<Real>(
24 : "relative_tolerance",
25 76 : 1.0e-3,
26 : "Relative tolerance on both compliance and cost to end the bisection method.");
27 76 : params.addRequiredCoupledVar("design_density", "Design density variable name.");
28 76 : params.addRequiredParam<VariableName>("density_sensitivity",
29 : "Name of the density_sensitivity variable.");
30 76 : params.addRequiredParam<VariableName>("cost_density_sensitivity",
31 : "Name of the cost density sensitivity variable.");
32 76 : params.addParam<VariableName>("thermal_sensitivity",
33 : "Name of the thermal density sensitivity variable.");
34 76 : params.addRequiredParam<Real>("volume_fraction", "Volume fraction.");
35 76 : params.addRequiredParam<Real>("cost_fraction", "Cost fraction.");
36 76 : params.addParam<Real>("bisection_move", 0.01, "Bisection move for the updated solution.");
37 76 : params.addParam<bool>("adaptive_move",
38 76 : false,
39 : "Whether incremental moves in the bisection algorithm will be reduced as "
40 : "the number of iterations increases. Note that the time must correspond to "
41 : "iteration number for better results.");
42 76 : params.addRequiredParam<VariableName>("cost", "Name of the cost variable.");
43 :
44 76 : params.addParam<Real>("bisection_lower_bound", 0, "Lower bound for the bisection algorithm.");
45 76 : params.addParam<Real>("bisection_upper_bound", 1e16, "Upper bound for the bisection algorithm.");
46 :
47 76 : params.addParam<std::vector<Real>>(
48 : "weight_mechanical_thermal",
49 : "List of values between 0 and 1 to weight the stiffness and thermal sensitivities");
50 :
51 76 : params.addParam<bool>("use_thermal_compliance",
52 76 : false,
53 : "Whether to include the thermal compliance in the sensitivities to "
54 : "minimize in conjunction with stiffness compliance.");
55 :
56 38 : return params;
57 0 : }
58 :
59 20 : DensityUpdateTwoConstraints::DensityUpdateTwoConstraints(const InputParameters & parameters)
60 : : ElementUserObject(parameters),
61 20 : _mesh(_subproblem.mesh()),
62 20 : _density_sensitivity_name(getParam<VariableName>("density_sensitivity")),
63 20 : _cost_density_sensitivity_name(getParam<VariableName>("cost_density_sensitivity")),
64 20 : _cost_name(getParam<VariableName>("cost")),
65 20 : _design_density(&writableVariable("design_density")),
66 20 : _density_sensitivity(&_subproblem.getStandardVariable(_tid, _density_sensitivity_name)),
67 20 : _cost_density_sensitivity(
68 20 : &_subproblem.getStandardVariable(_tid, _cost_density_sensitivity_name)),
69 20 : _cost(_subproblem.getStandardVariable(_tid, _cost_name)),
70 40 : _volume_fraction(getParam<Real>("volume_fraction")),
71 40 : _cost_fraction(getParam<Real>("cost_fraction")),
72 40 : _relative_tolerance(getParam<Real>("relative_tolerance")),
73 40 : _lower_bound(getParam<Real>("bisection_lower_bound")),
74 40 : _upper_bound(getParam<Real>("bisection_upper_bound")),
75 40 : _bisection_move(getParam<Real>("bisection_move")),
76 40 : _adaptive_move(getParam<bool>("adaptive_move")),
77 20 : _thermal_sensitivity_name(
78 40 : isParamValid("thermal_sensitivity") ? getParam<VariableName>("thermal_sensitivity") : ""),
79 40 : _thermal_sensitivity(isParamValid("thermal_sensitivity")
80 20 : ? &_subproblem.getStandardVariable(_tid, _thermal_sensitivity_name)
81 20 : : nullptr)
82 : {
83 40 : if (isParamValid("thermal_sensitivity"))
84 : {
85 0 : if (!isParamValid("thermal_sensitivity"))
86 0 : paramError("thermal_sensitivity",
87 : "The thermal_sensitivity variable name must be provided by the user if "
88 : "include_thermal_compliance is chosen to be true.");
89 :
90 0 : if (isParamValid("weight_mechanical_thermal"))
91 : {
92 0 : _weight_values = getParam<std::vector<Real>>("weight_mechanical_thermal");
93 0 : if (_weight_values.size() != 2)
94 0 : paramError("weight_mechanical_thermal",
95 : "Weighing of sensitivities is only available for the mechanical compliance and "
96 : "the thermal compliances problems, respectively.");
97 : }
98 : else
99 0 : paramError("weight_mechanical_thermal",
100 : "This parameter needs to be provided when including thermal sensitivity.");
101 : }
102 :
103 20 : auto transient = dynamic_cast<TransientBase *>(_app.getExecutioner());
104 :
105 20 : if (!transient && _adaptive_move)
106 0 : paramError("adaptive_move", "Cannot find a transient executioner for adaptive bisection move.");
107 :
108 20 : if (!dynamic_cast<MooseVariableFE<Real> *>(_design_density))
109 0 : paramError("design_density", "Design density must be a finite element variable");
110 :
111 20 : if (!dynamic_cast<const MooseVariableFE<Real> *>(_density_sensitivity))
112 0 : paramError("density_sensitivity", "Design sensitivity must be a finite element variable");
113 :
114 20 : if (!dynamic_cast<const MooseVariableFE<Real> *>(_cost_density_sensitivity))
115 0 : paramError("cost_density_sensitivity",
116 : "Cost density sensitivity must be a finite element variable");
117 20 : }
118 :
119 : void
120 56 : DensityUpdateTwoConstraints::timestepSetup()
121 : {
122 56 : gatherElementData();
123 56 : performOptimCritLoop();
124 56 : }
125 :
126 : void
127 11704 : DensityUpdateTwoConstraints::execute()
128 : {
129 : // Grab the element data for each id
130 11704 : auto elem_data_iter = _elem_data_map.find(_current_elem->id());
131 :
132 : // Check if the element data is not null
133 11704 : if (elem_data_iter != _elem_data_map.end())
134 : {
135 : ElementData & elem_data = elem_data_iter->second;
136 11704 : _design_density->setNodalValue(elem_data.new_density);
137 : }
138 : else
139 : {
140 0 : mooseError("Element data not found for the current element id.");
141 : }
142 11704 : }
143 :
144 : void
145 56 : DensityUpdateTwoConstraints::gatherElementData()
146 : {
147 : // Create parallel-consistent data structures constraining compliance and cost sensitivities.
148 : _elem_data_map.clear();
149 56 : _total_allowable_volume = 0;
150 56 : _total_allowable_cost = 0;
151 :
152 112 : for (const auto & sub_id : blockIDs())
153 13664 : for (const auto & elem : _mesh.getMesh().active_local_subdomain_elements_ptr_range(sub_id))
154 : {
155 13552 : dof_id_type elem_id = elem->id();
156 : ElementData data = ElementData(
157 13552 : dynamic_cast<MooseVariableFE<Real> *>(_design_density)->getElementalValue(elem),
158 13552 : dynamic_cast<const MooseVariableFE<Real> *>(_density_sensitivity)
159 : ->getElementalValue(elem),
160 13552 : dynamic_cast<const MooseVariableFE<Real> *>(_cost_density_sensitivity)
161 : ->getElementalValue(elem),
162 13552 : isParamValid("thermal_sensitivity") ? _thermal_sensitivity->getElementalValue(elem) : 0.0,
163 13552 : _cost.getElementalValue(elem),
164 13552 : elem->volume(),
165 67760 : 0);
166 13552 : _elem_data_map[elem_id] = data;
167 13552 : _total_allowable_volume += elem->volume();
168 56 : }
169 :
170 56 : _communicator.sum(_total_allowable_volume);
171 56 : _communicator.sum(_total_allowable_cost);
172 :
173 56 : _total_allowable_cost = _cost_fraction * _total_allowable_volume;
174 56 : _total_allowable_volume *= _volume_fraction;
175 56 : }
176 :
177 : void
178 56 : DensityUpdateTwoConstraints::performOptimCritLoop()
179 : {
180 : // Initialize the lower and upper bounds for the bisection method
181 56 : Real l1 = _lower_bound;
182 56 : Real l2 = _upper_bound;
183 :
184 : // Initialize the lower and upper bounds for cost solution
185 : Real c1 = _lower_bound;
186 : Real c2 = _upper_bound;
187 :
188 : bool perform_loop = true;
189 : // Loop until the relative difference between l1 and l2 is less than a small tolerance
190 :
191 63224 : while (perform_loop)
192 : {
193 : // Compute the midpoint between l1 and l2
194 63168 : Real lmid = 0.5 * (l2 + l1);
195 :
196 : // Compute the midpoint between c1 and c2
197 63168 : Real cmid = 0.5 * (c2 + c1);
198 :
199 : // Initialize the current total volume
200 63168 : Real curr_total_volume = 0;
201 :
202 : // Initialize the current total cost
203 63168 : Real curr_total_cost = 0;
204 :
205 : // Loop over all elements
206 15349824 : for (auto && [id, elem_data] : _elem_data_map)
207 : {
208 : // Compute the updated density for the current element
209 15286656 : Real new_density = computeUpdatedDensity(
210 : elem_data.old_density,
211 45859968 : elem_data.sensitivity * (isParamValid("thermal_sensitivity") ? _weight_values[0] : 1.0),
212 : elem_data.cost_sensitivity,
213 15286656 : elem_data.thermal_sensitivity *
214 30573312 : (isParamValid("thermal_sensitivity") ? _weight_values[1] : 1.0),
215 : elem_data.cost,
216 : lmid,
217 : cmid);
218 :
219 : // Update the current filtered density for the current element
220 15286656 : elem_data.new_density = new_density;
221 : // Update the current total volume
222 15286656 : curr_total_volume += new_density * elem_data.volume;
223 : // Update the current total cost
224 15286656 : curr_total_cost += new_density * elem_data.volume * elem_data.cost;
225 : }
226 :
227 : // Sum the current total volume across all processors
228 63168 : _communicator.sum(curr_total_volume);
229 63168 : _communicator.sum(curr_total_cost);
230 :
231 : // Update l1 or l2 based on whether the current total volume is greater than the total allowable
232 : // volume
233 63168 : if (curr_total_volume > _total_allowable_volume)
234 : l1 = lmid;
235 : else
236 : l2 = lmid;
237 :
238 63168 : if (curr_total_cost > _total_allowable_cost)
239 : c1 = cmid;
240 : else
241 : c2 = cmid;
242 :
243 : // Determine whether to continue the loop based on the relative difference between l1 and l2
244 : perform_loop =
245 63168 : (l2 - l1) / (l1 + l2) > _relative_tolerance || (c2 - c1) / (c1 + c2) > _relative_tolerance;
246 : }
247 56 : }
248 :
249 : // Method to compute the updated density for an element
250 : Real
251 15286656 : DensityUpdateTwoConstraints::computeUpdatedDensity(Real current_density,
252 : Real dc,
253 : Real cost_sensitivity,
254 : Real temp_sensitivity,
255 : Real cost,
256 : Real lmid,
257 : Real cmid)
258 : {
259 15286656 : Real move = _bisection_move;
260 :
261 : // Minimum move
262 15286656 : const Real move_min = 1.0e-3;
263 : // Control adaptivity
264 : const Real alpha = 0.96;
265 :
266 : // Adaptive move (takes the transient time step as SIMP iteration number)
267 15286656 : if (_adaptive_move)
268 0 : move = std::max(MathUtils::pow(alpha, _t) * move, move_min);
269 :
270 15286656 : Real denominator = lmid + cmid * cost + cmid * current_density * cost_sensitivity;
271 :
272 : // Effect of damping to be assessed
273 : const Real damping = 1.0;
274 :
275 : Real updated_density = std::max(
276 15286656 : 0.0,
277 : std::max(
278 15286656 : current_density - move,
279 15286656 : std::min(1.0,
280 30573312 : std::min(current_density + move,
281 30573312 : current_density *
282 15286656 : MathUtils::pow(std::sqrt((-(dc + temp_sensitivity) / denominator)),
283 15286656 : damping)))));
284 :
285 : // Return the updated density
286 15286656 : return updated_density;
287 : }
|
