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 424 : PenaltyWeightedGapUserObject::validParams()
21 : {
22 424 : InputParameters params = WeightedGapUserObject::validParams();
23 424 : params.addClassDescription("Computes the mortar normal contact force via a penalty approach.");
24 848 : params.addRequiredParam<Real>("penalty", "The penalty factor");
25 1272 : params.addRangeCheckedParam<Real>(
26 : "penalty_multiplier",
27 848 : 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 848 : params.addParam<bool>(
33 : "use_physical_gap",
34 848 : 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 848 : params.addRangeCheckedParam<Real>(
39 : "penetration_tolerance",
40 : "penetration_tolerance > 0",
41 : "Acceptable penetration distance at which augmented Lagrange iterations can be stopped");
42 1272 : params.addRangeCheckedParam<Real>(
43 : "max_penalty_multiplier",
44 848 : 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 848 : MooseEnum adaptivity_penalty_normal("SIMPLE BUSSETTA", "SIMPLE");
49 848 : adaptivity_penalty_normal.addDocumentation(
50 : "SIMPLE", "Keep multiplying by the penalty multiplier between AL iterations");
51 848 : 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 848 : params.addParam<MooseEnum>(
56 : "adaptivity_penalty_normal",
57 : adaptivity_penalty_normal,
58 : "The augmented Lagrange update strategy used on the normal penalty coefficient.");
59 848 : 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 848 : params.addParamNamesToGroup("penalty_multiplier penetration_tolerance max_penalty_multiplier",
64 : "Augmented Lagrange");
65 :
66 424 : return params;
67 424 : }
68 :
69 212 : PenaltyWeightedGapUserObject::PenaltyWeightedGapUserObject(const InputParameters & parameters)
70 : : WeightedGapUserObject(parameters),
71 : AugmentedLagrangeInterface(this),
72 212 : _penalty(getParam<Real>("penalty")),
73 424 : _penalty_multiplier(getParam<Real>("penalty_multiplier")),
74 212 : _penetration_tolerance(
75 450 : isParamValid("penetration_tolerance") ? getParam<Real>("penetration_tolerance") : 0.0),
76 212 : _augmented_lagrange_problem(
77 212 : dynamic_cast<AugmentedLagrangianContactProblemInterface *>(&_fe_problem)),
78 212 : _lagrangian_iteration_number(_augmented_lagrange_problem
79 212 : ? _augmented_lagrange_problem->getLagrangianIterationNumber()
80 : : _no_iterations),
81 212 : _dt(_fe_problem.dt()),
82 424 : _use_physical_gap(getParam<bool>("use_physical_gap")),
83 232 : _aux_lm_var(isCoupled("aux_lm") ? getVar("aux_lm", 0) : nullptr),
84 424 : _max_penalty_multiplier(getParam<Real>("max_penalty_multiplier")),
85 212 : _adaptivity_normal(
86 424 : getParam<MooseEnum>("adaptivity_penalty_normal").getEnum<AdaptivityNormalPenalty>())
87 : {
88 437 : auto check_type = [this](const auto & var, const auto & var_name)
89 : {
90 437 : 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 649 : };
95 212 : check_type(*_disp_x_var, "disp_x");
96 212 : check_type(*_disp_y_var, "disp_y");
97 212 : if (_has_disp_z)
98 13 : check_type(*_disp_z_var, "disp_z");
99 :
100 424 : 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 212 : }
105 :
106 : const VariableTestValue &
107 95 : PenaltyWeightedGapUserObject::test() const
108 : {
109 95 : return _aux_lm_var ? _aux_lm_var->phiLower() : _disp_x_var->phiLower();
110 : }
111 :
112 : const ADVariableValue &
113 7949488 : PenaltyWeightedGapUserObject::contactPressure() const
114 : {
115 7949488 : return _contact_pressure;
116 : }
117 :
118 : void
119 23715 : PenaltyWeightedGapUserObject::selfInitialize()
120 : {
121 117973 : for (auto & dof_pn : _dof_to_normal_pressure)
122 94258 : dof_pn.second = 0.0;
123 23715 : }
124 :
125 : void
126 8568 : PenaltyWeightedGapUserObject::initialize()
127 : {
128 8568 : WeightedGapUserObject::initialize();
129 8568 : selfInitialize();
130 8568 : }
131 :
132 : void
133 23715 : PenaltyWeightedGapUserObject::selfFinalize()
134 : {
135 : using std::min;
136 :
137 : // compute new normal pressure for each node
138 118945 : for (const auto & [dof_object, wgap] : _dof_to_weighted_gap)
139 : {
140 95230 : auto penalty = findValue(_dof_to_local_penalty, dof_object, _penalty);
141 :
142 : // If _use_physical_gap is true we "normalize" the penalty parameter with the surrounding area.
143 95230 : auto gap = _use_physical_gap ? adPhysicalGap(wgap) / wgap.second : wgap.first;
144 :
145 : const auto lagrange_multiplier =
146 144458 : _augmented_lagrange_problem ? _dof_to_lagrange_multiplier[dof_object] : 0.0;
147 :
148 : // keep the negative normal pressure (compressive per convention here)
149 95230 : auto normal_pressure = min(0.0, penalty * gap + lagrange_multiplier);
150 :
151 : // we switch conventins here and consider positive normal pressure as compressive
152 190460 : _dof_to_normal_pressure[dof_object] = -normal_pressure;
153 : }
154 23715 : }
155 :
156 : void
157 8568 : PenaltyWeightedGapUserObject::finalize()
158 : {
159 8568 : WeightedGapUserObject::finalize();
160 8568 : selfFinalize();
161 8568 : }
162 :
163 : Real
164 135304 : PenaltyWeightedGapUserObject::getNormalContactPressure(const Node * const node) const
165 : {
166 135304 : const auto it = _dof_to_normal_pressure.find(_subproblem.mesh().nodePtr(node->id()));
167 :
168 135304 : if (it != _dof_to_normal_pressure.end())
169 5716 : return MetaPhysicL::raw_value(it->second);
170 : else
171 : return 0.0;
172 : }
173 :
174 : Real
175 810 : PenaltyWeightedGapUserObject::getNormalLagrangeMultiplier(const Node * const node) const
176 : {
177 810 : const auto it = _dof_to_lagrange_multiplier.find(_subproblem.mesh().nodePtr(node->id()));
178 :
179 810 : if (it != _dof_to_lagrange_multiplier.end())
180 710 : return -MetaPhysicL::raw_value(it->second);
181 : else
182 : return 0.0;
183 : }
184 :
185 : void
186 169163 : PenaltyWeightedGapUserObject::reinit()
187 : {
188 169163 : _contact_pressure.resize(_qrule_msm->n_points());
189 638561 : for (const auto qp : make_range(_qrule_msm->n_points()))
190 469398 : _contact_pressure[qp] = 0.0;
191 :
192 638561 : for (const auto i : make_range(_test->size()))
193 : {
194 469398 : const Node * const node = _lower_secondary_elem->node_ptr(i);
195 :
196 1932482 : for (const auto qp : make_range(_qrule_msm->n_points()))
197 2926168 : _contact_pressure[qp] += (*_test)[i][qp] * _dof_to_normal_pressure[node];
198 : }
199 169163 : }
200 :
201 : void
202 1481 : PenaltyWeightedGapUserObject::selfTimestepSetup()
203 : {
204 : // Let's not clear the LMs for improved performance in
205 : // nonlinear problems
206 :
207 : // reset penalty
208 1552 : for (auto & dof_lp : _dof_to_local_penalty)
209 71 : dof_lp.second = _penalty;
210 :
211 : // clear previous gap
212 1552 : for (auto & dof_pg : _dof_to_previous_gap)
213 71 : dof_pg.second = 0.0;
214 :
215 : // save old timestep
216 1481 : _dt_old = _dt;
217 1481 : }
218 :
219 : void
220 681 : PenaltyWeightedGapUserObject::timestepSetup()
221 : {
222 681 : selfTimestepSetup();
223 681 : }
224 :
225 : bool
226 74 : PenaltyWeightedGapUserObject::isAugmentedLagrangianConverged()
227 : {
228 : using std::max, std::min, std::abs;
229 :
230 74 : Real max_positive_gap = 0.0;
231 74 : Real min_negative_gap = 0.0;
232 :
233 : // Get maximum gap to ascertain whether we are converged.
234 3788 : for (auto & [dof_object, wgap] : _dof_to_weighted_gap)
235 : {
236 3714 : const auto gap = physicalGap(wgap);
237 : {
238 : // Check condition for nodes that are active
239 7232 : if (gap < 0 || _dof_to_lagrange_multiplier[dof_object] < 0.0)
240 : {
241 282 : max_positive_gap = max(max_positive_gap, gap);
242 282 : min_negative_gap = min(min_negative_gap, gap);
243 : }
244 : }
245 : }
246 :
247 : // Communicate the extreme gap values in parallel.
248 : std::vector<Real> recv;
249 74 : if (this->_communicator.rank() == 0)
250 52 : recv.resize(this->_communicator.size());
251 74 : this->_communicator.gather(
252 74 : 0, -min_negative_gap > max_positive_gap ? min_negative_gap : max_positive_gap, recv);
253 :
254 74 : if (this->_communicator.rank() == 0)
255 : {
256 52 : min_negative_gap = *std::min_element(recv.begin(), recv.end());
257 52 : max_positive_gap = *std::max_element(recv.begin(), recv.end());
258 :
259 : // report the gap value with the largest magnitude
260 52 : const Real reported_gap =
261 52 : -min_negative_gap > max_positive_gap ? min_negative_gap : max_positive_gap;
262 52 : if (abs(reported_gap) > _penetration_tolerance)
263 : {
264 9 : mooseInfoRepeated("Penetration tolerance fail max_gap = ",
265 : reported_gap,
266 : " (gap_tol=",
267 9 : _penetration_tolerance,
268 : "). Iteration number is: ",
269 : _lagrangian_iteration_number,
270 : ".");
271 9 : return false;
272 : }
273 : else
274 43 : mooseInfoRepeated("Penetration tolerance success max_gap = ",
275 : reported_gap,
276 : " (gap_tol=",
277 43 : _penetration_tolerance,
278 : "). Iteration number is: ",
279 : _lagrangian_iteration_number,
280 : ".");
281 : }
282 :
283 : return true;
284 74 : }
285 :
286 : void
287 60 : PenaltyWeightedGapUserObject::augmentedLagrangianSetup()
288 : {
289 : // Loop over all nodes for which a gap has been computed
290 3060 : for (auto & [dof_object, wgap] : _dof_to_weighted_gap)
291 : {
292 : const Real gap = physicalGap(wgap);
293 : // Store previous augmented lagrange iteration gap
294 3000 : _dof_to_previous_gap[dof_object] = gap;
295 : }
296 60 : }
297 :
298 : void
299 60 : PenaltyWeightedGapUserObject::updateAugmentedLagrangianMultipliers()
300 : {
301 : using std::abs;
302 :
303 3060 : for (const auto & [dof_object, wgap] : _dof_to_weighted_gap)
304 : {
305 : auto & penalty = _dof_to_local_penalty[dof_object];
306 3000 : if (penalty == 0.0)
307 643 : penalty = _penalty;
308 :
309 3000 : const auto gap = getNormalGap(static_cast<const Node *>(dof_object));
310 : auto & lagrange_multiplier = _dof_to_lagrange_multiplier[dof_object];
311 :
312 : const auto possible_normalization =
313 3000 : (_use_physical_gap ? libmesh_map_find(_dof_to_weighted_gap, dof_object).second : 1.0);
314 :
315 : // Update penalty (the factor of 1/4 is suggested in the literature, the limit on AL iteration
316 : // caps the penalty increase)
317 : // Before we were updating the LM before adapting the penalty factor
318 3000 : if (lagrange_multiplier + gap * penalty / possible_normalization <= 0)
319 228 : lagrange_multiplier += gap * penalty / possible_normalization;
320 : else
321 2772 : lagrange_multiplier = 0.0;
322 :
323 3000 : const auto previous_gap = _dof_to_previous_gap[dof_object];
324 3000 : Real eval_tn = 0;
325 :
326 3000 : if (_adaptivity_normal == AdaptivityNormalPenalty::SIMPLE)
327 : {
328 0 : if (abs(gap) > 0.25 * abs(previous_gap) && abs(gap) > _penetration_tolerance)
329 0 : penalty *= _penalty_multiplier;
330 : }
331 3000 : else if (_adaptivity_normal == AdaptivityNormalPenalty::BUSSETTA)
332 3000 : bussettaAdaptivePenalty(previous_gap, gap, penalty, eval_tn);
333 :
334 : // Heuristics to bound the penalty factor
335 3000 : if (penalty < _penalty)
336 0 : penalty = _penalty;
337 3000 : else if (penalty > _max_penalty_multiplier * _penalty)
338 0 : penalty = _max_penalty_multiplier * _penalty;
339 : }
340 60 : }
341 :
342 : void
343 3000 : PenaltyWeightedGapUserObject::bussettaAdaptivePenalty(const Real previous_gap,
344 : const Real gap,
345 : Real & penalty,
346 : Real & eval_tn)
347 : {
348 : // Positive gaps means no contact
349 3000 : if (previous_gap > 0.0)
350 : {
351 2226 : penalty = _penalty;
352 2226 : eval_tn = 0.0;
353 : }
354 : else
355 : {
356 774 : if (previous_gap * gap < 0)
357 18 : eval_tn = 0.0;
358 : else
359 756 : eval_tn = penalty * previous_gap;
360 :
361 774 : adaptiveNormalPenalty(previous_gap, gap, penalty);
362 : }
363 3000 : }
364 :
365 : void
366 774 : PenaltyWeightedGapUserObject::adaptiveNormalPenalty(const Real previous_gap,
367 : const Real gap,
368 : Real & penalty)
369 : {
370 : using std::abs, std::max, std::sqrt;
371 :
372 774 : const bool condition_one = abs(abs(previous_gap / gap) - 1.0) < 0.01;
373 774 : const bool condition_two = abs(gap) > 1.01 * abs(previous_gap);
374 :
375 774 : if (previous_gap * gap < 0)
376 : {
377 18 : if (previous_gap > _penetration_tolerance)
378 0 : penalty = abs(penalty * previous_gap / gap * (abs(gap) + _penetration_tolerance) /
379 0 : (gap - previous_gap));
380 : else
381 18 : penalty = abs(penalty * previous_gap / 10.0 / gap);
382 : }
383 756 : else if (abs(gap) > _penetration_tolerance)
384 : {
385 643 : if (abs(gap - previous_gap) >
386 1335 : max(gap / 10.0, max(previous_gap / 10.0, 5.0 * _penetration_tolerance)))
387 643 : penalty *= 10.0;
388 0 : else if (condition_one && gap < 10.0 * _penetration_tolerance)
389 0 : penalty *= MathUtils::pow(sqrt(abs(gap) / _penetration_tolerance - 1.0) + 1.0, 2);
390 0 : else if (condition_two)
391 0 : penalty *= 10.0 * previous_gap / gap;
392 : else
393 0 : penalty *= sqrt(abs(gap) / _penetration_tolerance - 1.0) + 1.0;
394 : }
395 774 : }
|
