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 "PenaltyWeightedGapUserObject.h"
11 : #include "AugmentedLagrangianContactProblem.h"
12 : #include "MooseVariableFE.h"
13 : #include "SystemBase.h"
14 :
15 : registerMooseObject("ContactApp", PenaltyWeightedGapUserObject);
16 :
17 : const unsigned int PenaltyWeightedGapUserObject::_no_iterations = 0;
18 :
19 : InputParameters
20 684 : PenaltyWeightedGapUserObject::validParams()
21 : {
22 684 : InputParameters params = WeightedGapUserObject::validParams();
23 684 : params.addClassDescription("Computes the mortar normal contact force via a penalty approach.");
24 1368 : params.addRequiredParam<Real>("penalty", "The penalty factor");
25 2052 : params.addRangeCheckedParam<Real>(
26 : "penalty_multiplier",
27 1368 : 1.0,
28 : "penalty_multiplier > 0",
29 : "The penalty growth factor between augmented Lagrange "
30 : "iterations if weighted gap does not get closed fast enough. For frictional simulations, "
31 : "values smaller than 100 are recommended, e.g. 5.");
32 1368 : params.addParam<bool>(
33 : "use_physical_gap",
34 1368 : false,
35 : "Whether to use the physical normal gap (not scaled by mortar integration) and normalize the "
36 : "penalty coefficient with a representative lower-dimensional volume assigned to the node in "
37 : "the contacting boundary. This parameter is defaulted to 'true' in the contact action.");
38 1368 : params.addRangeCheckedParam<Real>(
39 : "penetration_tolerance",
40 : "penetration_tolerance > 0",
41 : "Acceptable penetration distance at which augmented Lagrange iterations can be stopped");
42 2052 : params.addRangeCheckedParam<Real>(
43 : "max_penalty_multiplier",
44 1368 : 1.0e3,
45 : "max_penalty_multiplier >= 1.0",
46 : "Maximum multiplier applied to penalty factors when adaptivity is used in an augmented "
47 : "Lagrange setting. The penalty factor supplied by the user is used as a reference.");
48 1368 : MooseEnum adaptivity_penalty_normal("SIMPLE BUSSETTA", "SIMPLE");
49 1368 : adaptivity_penalty_normal.addDocumentation(
50 : "SIMPLE", "Keep multiplying by the penalty multiplier between AL iterations");
51 1368 : adaptivity_penalty_normal.addDocumentation(
52 : "BUSSETTA",
53 : "Modify the penalty using an algorithm from Bussetta et al, 2012, Comput Mech 49:259-275 "
54 : "between AL iterations.");
55 1368 : params.addParam<MooseEnum>(
56 : "adaptivity_penalty_normal",
57 : adaptivity_penalty_normal,
58 : "The augmented Lagrange update strategy used on the normal penalty coefficient.");
59 1368 : params.addCoupledVar(
60 : "aux_lm",
61 : "Auxiliary variable that is utilized together with the "
62 : "penalty approach to interpolate the resulting contact tractions using dual bases.");
63 1368 : params.addParamNamesToGroup("penalty_multiplier penetration_tolerance max_penalty_multiplier",
64 : "Augmented Lagrange");
65 :
66 684 : return params;
67 684 : }
68 :
69 342 : PenaltyWeightedGapUserObject::PenaltyWeightedGapUserObject(const InputParameters & parameters)
70 : : WeightedGapUserObject(parameters),
71 : AugmentedLagrangeInterface(this),
72 342 : _penalty(getParam<Real>("penalty")),
73 684 : _penalty_multiplier(getParam<Real>("penalty_multiplier")),
74 342 : _penetration_tolerance(
75 756 : isParamValid("penetration_tolerance") ? getParam<Real>("penetration_tolerance") : 0.0),
76 342 : _augmented_lagrange_problem(
77 342 : dynamic_cast<AugmentedLagrangianContactProblemInterface *>(&_fe_problem)),
78 342 : _lagrangian_iteration_number(_augmented_lagrange_problem
79 342 : ? _augmented_lagrange_problem->getLagrangianIterationNumber()
80 : : _no_iterations),
81 342 : _dt(_fe_problem.dt()),
82 684 : _use_physical_gap(getParam<bool>("use_physical_gap")),
83 388 : _aux_lm_var(isCoupled("aux_lm") ? getVar("aux_lm", 0) : nullptr),
84 684 : _max_penalty_multiplier(getParam<Real>("max_penalty_multiplier")),
85 342 : _adaptivity_normal(
86 684 : getParam<MooseEnum>("adaptivity_penalty_normal").getEnum<AdaptivityNormalPenalty>())
87 : {
88 700 : auto check_type = [this](const auto & var, const auto & var_name)
89 : {
90 700 : if (!var.isNodal())
91 0 : paramError(var_name,
92 : "The displacement variables must have degrees of freedom exclusively on "
93 : "nodes, e.g. they should probably be of finite element type 'Lagrange'.");
94 1042 : };
95 342 : check_type(*_disp_x_var, "disp_x");
96 342 : check_type(*_disp_y_var, "disp_y");
97 342 : if (_has_disp_z)
98 16 : check_type(*_disp_z_var, "disp_z");
99 :
100 684 : if (!_augmented_lagrange_problem == isParamValid("penetration_tolerance"))
101 0 : paramError("penetration_tolerance",
102 : "This parameter must be supplied if and only if an augmented Lagrange problem "
103 : "object is used.");
104 342 : }
105 :
106 : const VariableTestValue &
107 171 : PenaltyWeightedGapUserObject::test() const
108 : {
109 171 : return _aux_lm_var ? _aux_lm_var->phiLower() : _disp_x_var->phiLower();
110 : }
111 :
112 : const ADVariableValue &
113 8673472 : PenaltyWeightedGapUserObject::contactPressure() const
114 : {
115 8673472 : return _contact_pressure;
116 : }
117 :
118 : void
119 29775 : PenaltyWeightedGapUserObject::selfInitialize()
120 : {
121 162846 : for (auto & dof_pn : _dof_to_normal_pressure)
122 133071 : dof_pn.second = 0.0;
123 29775 : }
124 :
125 : void
126 13362 : PenaltyWeightedGapUserObject::initialize()
127 : {
128 13362 : WeightedGapUserObject::initialize();
129 13362 : selfInitialize();
130 13362 : }
131 :
132 : void
133 29775 : PenaltyWeightedGapUserObject::selfFinalize()
134 : {
135 : // compute new normal pressure for each node
136 164440 : for (const auto & [dof_object, wgap] : _dof_to_weighted_gap)
137 : {
138 134665 : auto penalty = findValue(_dof_to_local_penalty, dof_object, _penalty);
139 :
140 : // If _use_physical_gap is true we "normalize" the penalty parameter with the surrounding area.
141 134665 : auto gap = _use_physical_gap ? adPhysicalGap(wgap) / wgap.second : wgap.first;
142 :
143 : const auto lagrange_multiplier =
144 219373 : _augmented_lagrange_problem ? _dof_to_lagrange_multiplier[dof_object] : 0.0;
145 :
146 : // keep the negative normal pressure (compressive per convention here)
147 134665 : auto normal_pressure = std::min(0.0, penalty * gap + lagrange_multiplier);
148 :
149 : // we switch conventins here and consider positive normal pressure as compressive
150 269330 : _dof_to_normal_pressure[dof_object] = -normal_pressure;
151 : }
152 29775 : }
153 :
154 : void
155 13362 : PenaltyWeightedGapUserObject::finalize()
156 : {
157 13362 : WeightedGapUserObject::finalize();
158 13362 : selfFinalize();
159 13362 : }
160 :
161 : Real
162 129696 : PenaltyWeightedGapUserObject::getNormalContactPressure(const Node * const node) const
163 : {
164 129696 : const auto it = _dof_to_normal_pressure.find(_subproblem.mesh().nodePtr(node->id()));
165 :
166 129696 : if (it != _dof_to_normal_pressure.end())
167 4136 : return MetaPhysicL::raw_value(it->second);
168 : else
169 : return 0.0;
170 : }
171 :
172 : Real
173 1296 : PenaltyWeightedGapUserObject::getNormalLagrangeMultiplier(const Node * const node) const
174 : {
175 1296 : const auto it = _dof_to_lagrange_multiplier.find(_subproblem.mesh().nodePtr(node->id()));
176 :
177 1296 : if (it != _dof_to_lagrange_multiplier.end())
178 1136 : return -MetaPhysicL::raw_value(it->second);
179 : else
180 : return 0.0;
181 : }
182 :
183 : void
184 209292 : PenaltyWeightedGapUserObject::reinit()
185 : {
186 209292 : _contact_pressure.resize(_qrule_msm->n_points());
187 760996 : for (const auto qp : make_range(_qrule_msm->n_points()))
188 551704 : _contact_pressure[qp] = 0.0;
189 :
190 760996 : for (const auto i : make_range(_test->size()))
191 : {
192 551704 : const Node * const node = _lower_secondary_elem->node_ptr(i);
193 :
194 2187592 : for (const auto qp : make_range(_qrule_msm->n_points()))
195 3271776 : _contact_pressure[qp] += (*_test)[i][qp] * _dof_to_normal_pressure[node];
196 : }
197 209292 : }
198 :
199 : void
200 1990 : PenaltyWeightedGapUserObject::selfTimestepSetup()
201 : {
202 : // Let's not clear the LMs for improved performance in
203 : // nonlinear problems
204 :
205 : // reset penalty
206 1990 : for (auto & dof_lp : _dof_to_local_penalty)
207 0 : dof_lp.second = _penalty;
208 :
209 : // clear previous gap
210 1990 : for (auto & dof_pg : _dof_to_previous_gap)
211 0 : dof_pg.second = 0.0;
212 :
213 : // save old timestep
214 1990 : _dt_old = _dt;
215 1990 : }
216 :
217 : void
218 1062 : PenaltyWeightedGapUserObject::timestepSetup()
219 : {
220 1062 : selfTimestepSetup();
221 1062 : }
222 :
223 : bool
224 143 : PenaltyWeightedGapUserObject::isAugmentedLagrangianConverged()
225 : {
226 143 : Real max_positive_gap = 0.0;
227 143 : Real min_negative_gap = 0.0;
228 :
229 : // Get maximum gap to ascertain whether we are converged.
230 6446 : for (auto & [dof_object, wgap] : _dof_to_weighted_gap)
231 : {
232 6303 : const auto gap = physicalGap(wgap);
233 : {
234 : // Check condition for nodes that are active
235 12255 : if (gap < 0 || _dof_to_lagrange_multiplier[dof_object] < 0.0)
236 : {
237 495 : max_positive_gap = std::max(max_positive_gap, gap);
238 495 : min_negative_gap = std::min(min_negative_gap, gap);
239 : }
240 : }
241 : }
242 :
243 : // Communicate the extreme gap values in parallel.
244 : std::vector<Real> recv;
245 143 : if (this->_communicator.rank() == 0)
246 88 : recv.resize(this->_communicator.size());
247 143 : this->_communicator.gather(
248 143 : 0, -min_negative_gap > max_positive_gap ? min_negative_gap : max_positive_gap, recv);
249 :
250 143 : if (this->_communicator.rank() == 0)
251 : {
252 88 : min_negative_gap = *std::min_element(recv.begin(), recv.end());
253 88 : max_positive_gap = *std::max_element(recv.begin(), recv.end());
254 :
255 : // report the gap value with the largest magnitude
256 88 : const Real reported_gap =
257 88 : -min_negative_gap > max_positive_gap ? min_negative_gap : max_positive_gap;
258 88 : if (std::abs(reported_gap) > _penetration_tolerance)
259 : {
260 16 : mooseInfoRepeated("Penetration tolerance fail max_gap = ",
261 : reported_gap,
262 : " (gap_tol=",
263 16 : _penetration_tolerance,
264 : "). Iteration number is: ",
265 : _lagrangian_iteration_number,
266 : ".");
267 16 : return false;
268 : }
269 : else
270 72 : mooseInfoRepeated("Penetration tolerance success max_gap = ",
271 : reported_gap,
272 : " (gap_tol=",
273 72 : _penetration_tolerance,
274 : "). Iteration number is: ",
275 : _lagrangian_iteration_number,
276 : ".");
277 : }
278 :
279 : return true;
280 : }
281 :
282 : void
283 117 : PenaltyWeightedGapUserObject::augmentedLagrangianSetup()
284 : {
285 : // Loop over all nodes for which a gap has been computed
286 5274 : for (auto & [dof_object, wgap] : _dof_to_weighted_gap)
287 : {
288 : const Real gap = physicalGap(wgap);
289 : // Store previous augmented lagrange iteration gap
290 5157 : _dof_to_previous_gap[dof_object] = gap;
291 : }
292 117 : }
293 :
294 : void
295 117 : PenaltyWeightedGapUserObject::updateAugmentedLagrangianMultipliers()
296 : {
297 5274 : for (const auto & [dof_object, wgap] : _dof_to_weighted_gap)
298 : {
299 : auto & penalty = _dof_to_local_penalty[dof_object];
300 5157 : if (penalty == 0.0)
301 1146 : penalty = _penalty;
302 :
303 5157 : const auto gap = getNormalGap(static_cast<const Node *>(dof_object));
304 : auto & lagrange_multiplier = _dof_to_lagrange_multiplier[dof_object];
305 :
306 : const auto possible_normalization =
307 5157 : (_use_physical_gap ? libmesh_map_find(_dof_to_weighted_gap, dof_object).second : 1.0);
308 :
309 : // Update penalty (the factor of 1/4 is suggested in the literature, the limit on AL iteration
310 : // caps the penalty increase)
311 : // Before we were updating the LM before adapting the penalty factor
312 5157 : if (lagrange_multiplier + gap * penalty / possible_normalization <= 0)
313 405 : lagrange_multiplier += gap * penalty / possible_normalization;
314 : else
315 4752 : lagrange_multiplier = 0.0;
316 :
317 5157 : const auto previous_gap = _dof_to_previous_gap[dof_object];
318 5157 : Real eval_tn = 0;
319 :
320 5157 : if (_adaptivity_normal == AdaptivityNormalPenalty::SIMPLE)
321 : {
322 0 : if (std::abs(gap) > 0.25 * std::abs(previous_gap) && std::abs(gap) > _penetration_tolerance)
323 0 : penalty *= _penalty_multiplier;
324 : }
325 5157 : else if (_adaptivity_normal == AdaptivityNormalPenalty::BUSSETTA)
326 5157 : bussettaAdaptivePenalty(previous_gap, gap, penalty, eval_tn);
327 :
328 : // Heuristics to bound the penalty factor
329 5157 : if (penalty < _penalty)
330 0 : penalty = _penalty;
331 5157 : else if (penalty > _max_penalty_multiplier * _penalty)
332 0 : penalty = _max_penalty_multiplier * _penalty;
333 : }
334 117 : }
335 :
336 : void
337 5157 : PenaltyWeightedGapUserObject::bussettaAdaptivePenalty(const Real previous_gap,
338 : const Real gap,
339 : Real & penalty,
340 : Real & eval_tn)
341 : {
342 : // Positive gaps means no contact
343 5157 : if (previous_gap > 0.0)
344 : {
345 3776 : penalty = _penalty;
346 3776 : eval_tn = 0.0;
347 : }
348 : else
349 : {
350 1381 : if (previous_gap * gap < 0)
351 32 : eval_tn = 0.0;
352 : else
353 1349 : eval_tn = penalty * previous_gap;
354 :
355 1381 : adaptiveNormalPenalty(previous_gap, gap, penalty);
356 : }
357 5157 : }
358 :
359 : void
360 1381 : PenaltyWeightedGapUserObject::adaptiveNormalPenalty(const Real previous_gap,
361 : const Real gap,
362 : Real & penalty)
363 : {
364 1381 : const bool condition_one = std::abs(std::abs(previous_gap / gap) - 1.0) < 0.01;
365 1381 : const bool condition_two = std::abs(gap) > 1.01 * std::abs(previous_gap);
366 :
367 1381 : if (previous_gap * gap < 0)
368 : {
369 32 : if (previous_gap > _penetration_tolerance)
370 0 : penalty = std::abs(penalty * previous_gap / gap * (std::abs(gap) + _penetration_tolerance) /
371 0 : (gap - previous_gap));
372 : else
373 32 : penalty = std::abs(penalty * previous_gap / 10.0 / gap);
374 : }
375 1349 : else if (std::abs(gap) > _penetration_tolerance)
376 : {
377 1146 : if (std::abs(gap - previous_gap) >
378 2382 : std::max(gap / 10.0, std::max(previous_gap / 10.0, 5.0 * _penetration_tolerance)))
379 1146 : penalty *= 10.0;
380 0 : else if (condition_one && gap < 10.0 * _penetration_tolerance)
381 0 : penalty *= MathUtils::pow(std::sqrt(std::abs(gap) / _penetration_tolerance - 1.0) + 1.0, 2);
382 0 : else if (condition_two)
383 0 : penalty *= 10.0 * previous_gap / gap;
384 : else
385 0 : penalty *= std::sqrt(std::abs(gap) / _penetration_tolerance - 1.0) + 1.0;
386 : }
387 1381 : }
|
