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 82 : DensityUpdateTwoConstraints::validParams()
18 : {
19 82 : InputParameters params = ElementUserObject::validParams();
20 82 : params.addClassDescription(
21 : "Compute updated densities based on sensitivities using an optimality criteria method to "
22 : "keep the volume and cost constraints satisified.");
23 164 : params.addParam<Real>(
24 : "relative_tolerance",
25 164 : 1.0e-3,
26 : "Relative tolerance on both compliance and cost to end the bisection method.");
27 164 : params.addRequiredCoupledVar("design_density", "Design density variable name.");
28 164 : params.addRequiredParam<VariableName>("density_sensitivity",
29 : "Name of the density_sensitivity variable.");
30 164 : params.addRequiredParam<VariableName>("cost_density_sensitivity",
31 : "Name of the cost density sensitivity variable.");
32 164 : params.addParam<VariableName>("thermal_sensitivity",
33 : "Name of the thermal density sensitivity variable.");
34 164 : params.addRequiredParam<Real>("volume_fraction", "Volume fraction.");
35 164 : params.addRequiredParam<Real>("cost_fraction", "Cost fraction.");
36 164 : params.addParam<Real>("bisection_move", 0.01, "Bisection move for the updated solution.");
37 164 : params.addParam<bool>("adaptive_move",
38 164 : 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 164 : params.addRequiredParam<VariableName>("cost", "Name of the cost variable.");
43 :
44 164 : params.addParam<Real>("bisection_lower_bound", 0, "Lower bound for the bisection algorithm.");
45 164 : params.addParam<Real>("bisection_upper_bound", 1e16, "Upper bound for the bisection algorithm.");
46 :
47 164 : 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 164 : params.addParam<bool>("use_thermal_compliance",
52 164 : false,
53 : "Whether to include the thermal compliance in the sensitivities to "
54 : "minimize in conjunction with stiffness compliance.");
55 :
56 82 : return params;
57 0 : }
58 :
59 44 : DensityUpdateTwoConstraints::DensityUpdateTwoConstraints(const InputParameters & parameters)
60 : : ElementUserObject(parameters),
61 44 : _mesh(_subproblem.mesh()),
62 44 : _density_sensitivity_name(getParam<VariableName>("density_sensitivity")),
63 44 : _cost_density_sensitivity_name(getParam<VariableName>("cost_density_sensitivity")),
64 44 : _cost_name(getParam<VariableName>("cost")),
65 44 : _design_density(&writableVariable("design_density")),
66 44 : _density_sensitivity(&_subproblem.getStandardVariable(_tid, _density_sensitivity_name)),
67 44 : _cost_density_sensitivity(
68 44 : &_subproblem.getStandardVariable(_tid, _cost_density_sensitivity_name)),
69 44 : _cost(_subproblem.getStandardVariable(_tid, _cost_name)),
70 88 : _volume_fraction(getParam<Real>("volume_fraction")),
71 88 : _cost_fraction(getParam<Real>("cost_fraction")),
72 88 : _relative_tolerance(getParam<Real>("relative_tolerance")),
73 88 : _lower_bound(getParam<Real>("bisection_lower_bound")),
74 88 : _upper_bound(getParam<Real>("bisection_upper_bound")),
75 88 : _bisection_move(getParam<Real>("bisection_move")),
76 88 : _adaptive_move(getParam<bool>("adaptive_move")),
77 44 : _thermal_sensitivity_name(
78 88 : isParamValid("thermal_sensitivity") ? getParam<VariableName>("thermal_sensitivity") : ""),
79 88 : _thermal_sensitivity(isParamValid("thermal_sensitivity")
80 44 : ? &_subproblem.getStandardVariable(_tid, _thermal_sensitivity_name)
81 44 : : nullptr)
82 : {
83 88 : 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 44 : auto transient = dynamic_cast<TransientBase *>(_app.getExecutioner());
104 :
105 44 : if (!transient && _adaptive_move)
106 0 : paramError("adaptive_move", "Cannot find a transient executioner for adaptive bisection move.");
107 :
108 44 : if (!dynamic_cast<MooseVariableFE<Real> *>(_design_density))
109 0 : paramError("design_density", "Design density must be a finite element variable");
110 :
111 44 : if (!dynamic_cast<const MooseVariableFE<Real> *>(_density_sensitivity))
112 0 : paramError("density_sensitivity", "Design sensitivity must be a finite element variable");
113 :
114 44 : 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 44 : }
118 :
119 : void
120 104 : DensityUpdateTwoConstraints::timestepSetup()
121 : {
122 104 : gatherElementData();
123 104 : performOptimCritLoop();
124 104 : }
125 :
126 : void
127 17248 : DensityUpdateTwoConstraints::execute()
128 : {
129 : // Grab the element data for each id
130 17248 : auto elem_data_iter = _elem_data_map.find(_current_elem->id());
131 :
132 : // Check if the element data is not null
133 17248 : if (elem_data_iter != _elem_data_map.end())
134 : {
135 : ElementData & elem_data = elem_data_iter->second;
136 17248 : _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 17248 : }
143 :
144 : void
145 104 : DensityUpdateTwoConstraints::gatherElementData()
146 : {
147 : // Create parallel-consistent data structures constraining compliance and cost sensitivities.
148 : _elem_data_map.clear();
149 104 : _total_allowable_volume = 0;
150 104 : _total_allowable_cost = 0;
151 :
152 208 : for (const auto & sub_id : blockIDs())
153 21152 : for (const auto & elem : _mesh.getMesh().active_local_subdomain_elements_ptr_range(sub_id))
154 : {
155 20944 : dof_id_type elem_id = elem->id();
156 : ElementData data = ElementData(
157 20944 : dynamic_cast<MooseVariableFE<Real> *>(_design_density)->getElementalValue(elem),
158 20944 : dynamic_cast<const MooseVariableFE<Real> *>(_density_sensitivity)
159 : ->getElementalValue(elem),
160 20944 : dynamic_cast<const MooseVariableFE<Real> *>(_cost_density_sensitivity)
161 : ->getElementalValue(elem),
162 20944 : isParamValid("thermal_sensitivity") ? _thermal_sensitivity->getElementalValue(elem) : 0.0,
163 20944 : _cost.getElementalValue(elem),
164 20944 : elem->volume(),
165 104720 : 0);
166 20944 : _elem_data_map[elem_id] = data;
167 20944 : _total_allowable_volume += elem->volume();
168 104 : }
169 :
170 104 : _communicator.sum(_total_allowable_volume);
171 104 : _communicator.sum(_total_allowable_cost);
172 :
173 104 : _total_allowable_cost = _cost_fraction * _total_allowable_volume;
174 104 : _total_allowable_volume *= _volume_fraction;
175 104 : }
176 :
177 : void
178 104 : DensityUpdateTwoConstraints::performOptimCritLoop()
179 : {
180 : // Initialize the lower and upper bounds for the bisection method
181 104 : Real l1 = _lower_bound;
182 104 : 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 117416 : while (perform_loop)
192 : {
193 : // Compute the midpoint between l1 and l2
194 117312 : Real lmid = 0.5 * (l2 + l1);
195 :
196 : // Compute the midpoint between c1 and c2
197 117312 : Real cmid = 0.5 * (c2 + c1);
198 :
199 : // Initialize the current total volume
200 117312 : Real curr_total_volume = 0;
201 :
202 : // Initialize the current total cost
203 117312 : Real curr_total_cost = 0;
204 :
205 : // Loop over all elements
206 23742144 : for (auto && [id, elem_data] : _elem_data_map)
207 : {
208 : // Compute the updated density for the current element
209 23624832 : Real new_density = computeUpdatedDensity(
210 : elem_data.old_density,
211 70874496 : elem_data.sensitivity * (isParamValid("thermal_sensitivity") ? _weight_values[0] : 1.0),
212 : elem_data.cost_sensitivity,
213 23624832 : elem_data.thermal_sensitivity *
214 47249664 : (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 23624832 : elem_data.new_density = new_density;
221 : // Update the current total volume
222 23624832 : curr_total_volume += new_density * elem_data.volume;
223 : // Update the current total cost
224 23624832 : curr_total_cost += new_density * elem_data.volume * elem_data.cost;
225 : }
226 :
227 : // Sum the current total volume across all processors
228 117312 : _communicator.sum(curr_total_volume);
229 117312 : _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 117312 : if (curr_total_volume > _total_allowable_volume)
234 : l1 = lmid;
235 : else
236 : l2 = lmid;
237 :
238 117312 : 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 117312 : (l2 - l1) / (l1 + l2) > _relative_tolerance || (c2 - c1) / (c1 + c2) > _relative_tolerance;
246 : }
247 104 : }
248 :
249 : // Method to compute the updated density for an element
250 : Real
251 23624832 : 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 23624832 : Real move = _bisection_move;
260 :
261 : // Minimum move
262 23624832 : 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 23624832 : if (_adaptive_move)
268 0 : move = std::max(MathUtils::pow(alpha, _t) * move, move_min);
269 :
270 23624832 : 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 23624832 : 0.0,
277 : std::max(
278 23624832 : current_density - move,
279 23624832 : std::min(1.0,
280 47249664 : std::min(current_density + move,
281 47249664 : current_density *
282 23624832 : MathUtils::pow(std::sqrt((-(dc + temp_sensitivity) / denominator)),
283 23624832 : damping)))));
284 :
285 : // Return the updated density
286 23624832 : return updated_density;
287 : }
|
