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 662 : PenaltyWeightedGapUserObject::validParams()
21 : {
22 662 : InputParameters params = WeightedGapUserObject::validParams();
23 662 : params.addClassDescription("Computes the mortar normal contact force via a penalty approach.");
24 1324 : params.addRequiredParam<Real>("penalty", "The penalty factor");
25 1986 : params.addRangeCheckedParam<Real>(
26 : "penalty_multiplier",
27 1324 : 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 1324 : params.addParam<bool>(
33 : "use_physical_gap",
34 1324 : 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 1324 : params.addRangeCheckedParam<Real>(
39 : "penetration_tolerance",
40 : "penetration_tolerance > 0",
41 : "Acceptable penetration distance at which augmented Lagrange iterations can be stopped");
42 1986 : params.addRangeCheckedParam<Real>(
43 : "max_penalty_multiplier",
44 1324 : 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 1324 : MooseEnum adaptivity_penalty_normal("SIMPLE BUSSETTA", "SIMPLE");
49 1324 : adaptivity_penalty_normal.addDocumentation(
50 : "SIMPLE", "Keep multiplying by the penalty multiplier between AL iterations");
51 1324 : 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 1324 : params.addParam<MooseEnum>(
56 : "adaptivity_penalty_normal",
57 : adaptivity_penalty_normal,
58 : "The augmented Lagrange update strategy used on the normal penalty coefficient.");
59 1324 : 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 1324 : params.addParamNamesToGroup("penalty_multiplier penetration_tolerance max_penalty_multiplier",
64 : "Augmented Lagrange");
65 :
66 662 : return params;
67 662 : }
68 :
69 331 : PenaltyWeightedGapUserObject::PenaltyWeightedGapUserObject(const InputParameters & parameters)
70 : : WeightedGapUserObject(parameters),
71 : AugmentedLagrangeInterface(this),
72 331 : _penalty(getParam<Real>("penalty")),
73 662 : _penalty_multiplier(getParam<Real>("penalty_multiplier")),
74 331 : _penetration_tolerance(
75 696 : isParamValid("penetration_tolerance") ? getParam<Real>("penetration_tolerance") : 0.0),
76 331 : _augmented_lagrange_problem(
77 331 : dynamic_cast<AugmentedLagrangianContactProblemInterface *>(&_fe_problem)),
78 331 : _lagrangian_iteration_number(_augmented_lagrange_problem
79 331 : ? _augmented_lagrange_problem->getLagrangianIterationNumber()
80 : : _no_iterations),
81 331 : _dt(_fe_problem.dt()),
82 662 : _use_physical_gap(getParam<bool>("use_physical_gap")),
83 359 : _aux_lm_var(isCoupled("aux_lm") ? getVar("aux_lm", 0) : nullptr),
84 662 : _max_penalty_multiplier(getParam<Real>("max_penalty_multiplier")),
85 331 : _adaptivity_normal(
86 662 : getParam<MooseEnum>("adaptivity_penalty_normal").getEnum<AdaptivityNormalPenalty>())
87 : {
88 679 : auto check_type = [this](const auto & var, const auto & var_name)
89 : {
90 679 : 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 1010 : };
95 331 : check_type(*_disp_x_var, "disp_x");
96 331 : check_type(*_disp_y_var, "disp_y");
97 331 : if (_has_disp_z)
98 17 : check_type(*_disp_z_var, "disp_z");
99 :
100 662 : 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 331 : }
105 :
106 : const VariableTestValue &
107 173 : PenaltyWeightedGapUserObject::test() const
108 : {
109 173 : return _aux_lm_var ? _aux_lm_var->phiLower() : _disp_x_var->phiLower();
110 : }
111 :
112 : const ADVariableValue &
113 9617920 : PenaltyWeightedGapUserObject::contactPressure() const
114 : {
115 9617920 : return _contact_pressure;
116 : }
117 :
118 : void
119 30120 : PenaltyWeightedGapUserObject::selfInitialize()
120 : {
121 145141 : for (auto & dof_pn : _dof_to_normal_pressure)
122 115021 : dof_pn.second = 0.0;
123 30120 : }
124 :
125 : void
126 13424 : PenaltyWeightedGapUserObject::initialize()
127 : {
128 13424 : WeightedGapUserObject::initialize();
129 13424 : selfInitialize();
130 13424 : }
131 :
132 : void
133 30120 : PenaltyWeightedGapUserObject::selfFinalize()
134 : {
135 : using std::min;
136 :
137 : // compute new normal pressure for each node
138 146402 : for (const auto & [dof_object, wgap] : _dof_to_weighted_gap)
139 : {
140 116282 : 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 116282 : auto gap = _use_physical_gap ? adPhysicalGap(wgap) / wgap.second : wgap.first;
144 :
145 : const auto lagrange_multiplier =
146 176154 : _augmented_lagrange_problem ? _dof_to_lagrange_multiplier[dof_object] : 0.0;
147 :
148 : // keep the negative normal pressure (compressive per convention here)
149 116282 : auto normal_pressure = min(0.0, penalty * gap + lagrange_multiplier);
150 :
151 : // we switch conventins here and consider positive normal pressure as compressive
152 232564 : _dof_to_normal_pressure[dof_object] = -normal_pressure;
153 : }
154 30120 : }
155 :
156 : void
157 13424 : PenaltyWeightedGapUserObject::finalize()
158 : {
159 13424 : WeightedGapUserObject::finalize();
160 13424 : selfFinalize();
161 13424 : }
162 :
163 : Real
164 160612 : PenaltyWeightedGapUserObject::getNormalContactPressure(const Node * const node) const
165 : {
166 160612 : const auto it = _dof_to_normal_pressure.find(_subproblem.mesh().nodePtr(node->id()));
167 :
168 160612 : if (it != _dof_to_normal_pressure.end())
169 6442 : return MetaPhysicL::raw_value(it->second);
170 : else
171 : return 0.0;
172 : }
173 :
174 : Real
175 972 : PenaltyWeightedGapUserObject::getNormalLagrangeMultiplier(const Node * const node) const
176 : {
177 972 : const auto it = _dof_to_lagrange_multiplier.find(_subproblem.mesh().nodePtr(node->id()));
178 :
179 972 : if (it != _dof_to_lagrange_multiplier.end())
180 852 : return -MetaPhysicL::raw_value(it->second);
181 : else
182 : return 0.0;
183 : }
184 :
185 : void
186 206880 : PenaltyWeightedGapUserObject::reinit()
187 : {
188 206880 : _contact_pressure.resize(_qrule_msm->n_points());
189 778336 : for (const auto qp : make_range(_qrule_msm->n_points()))
190 571456 : _contact_pressure[qp] = 0.0;
191 :
192 778336 : for (const auto i : make_range(_test->size()))
193 : {
194 571456 : const Node * const node = _lower_secondary_elem->node_ptr(i);
195 :
196 2345152 : for (const auto qp : make_range(_qrule_msm->n_points()))
197 3547392 : _contact_pressure[qp] += (*_test)[i][qp] * _dof_to_normal_pressure[node];
198 : }
199 206880 : }
200 :
201 : void
202 2062 : PenaltyWeightedGapUserObject::selfTimestepSetup()
203 : {
204 : // Let's not clear the LMs for improved performance in
205 : // nonlinear problems
206 :
207 : // reset penalty
208 2133 : for (auto & dof_lp : _dof_to_local_penalty)
209 71 : dof_lp.second = _penalty;
210 :
211 : // clear previous gap
212 2133 : for (auto & dof_pg : _dof_to_previous_gap)
213 71 : dof_pg.second = 0.0;
214 :
215 : // save old timestep
216 2062 : _dt_old = _dt;
217 2062 : }
218 :
219 : void
220 1067 : PenaltyWeightedGapUserObject::timestepSetup()
221 : {
222 1067 : selfTimestepSetup();
223 1067 : }
224 :
225 : bool
226 96 : PenaltyWeightedGapUserObject::isAugmentedLagrangianConverged()
227 : {
228 : using std::max, std::min, std::abs;
229 :
230 96 : Real max_positive_gap = 0.0;
231 96 : Real min_negative_gap = 0.0;
232 :
233 : // Get maximum gap to ascertain whether we are converged.
234 4602 : for (auto & [dof_object, wgap] : _dof_to_weighted_gap)
235 : {
236 4506 : const auto gap = physicalGap(wgap);
237 : {
238 : // Check condition for nodes that are active
239 8768 : if (gap < 0 || _dof_to_lagrange_multiplier[dof_object] < 0.0)
240 : {
241 348 : max_positive_gap = max(max_positive_gap, gap);
242 348 : 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 96 : if (this->_communicator.rank() == 0)
250 63 : recv.resize(this->_communicator.size());
251 96 : this->_communicator.gather(
252 96 : 0, -min_negative_gap > max_positive_gap ? min_negative_gap : max_positive_gap, recv);
253 :
254 96 : if (this->_communicator.rank() == 0)
255 : {
256 63 : min_negative_gap = *std::min_element(recv.begin(), recv.end());
257 63 : max_positive_gap = *std::max_element(recv.begin(), recv.end());
258 :
259 : // report the gap value with the largest magnitude
260 63 : const Real reported_gap =
261 63 : -min_negative_gap > max_positive_gap ? min_negative_gap : max_positive_gap;
262 63 : if (abs(reported_gap) > _penetration_tolerance)
263 : {
264 11 : mooseInfoRepeated("Penetration tolerance fail max_gap = ",
265 : reported_gap,
266 : " (gap_tol=",
267 11 : _penetration_tolerance,
268 : "). Iteration number is: ",
269 : _lagrangian_iteration_number,
270 : ".");
271 11 : return false;
272 : }
273 : else
274 52 : mooseInfoRepeated("Penetration tolerance success max_gap = ",
275 : reported_gap,
276 : " (gap_tol=",
277 52 : _penetration_tolerance,
278 : "). Iteration number is: ",
279 : _lagrangian_iteration_number,
280 : ".");
281 : }
282 :
283 : return true;
284 96 : }
285 :
286 : void
287 78 : PenaltyWeightedGapUserObject::augmentedLagrangianSetup()
288 : {
289 : // Loop over all nodes for which a gap has been computed
290 3726 : for (auto & [dof_object, wgap] : _dof_to_weighted_gap)
291 : {
292 : const Real gap = physicalGap(wgap);
293 : // Store previous augmented lagrange iteration gap
294 3648 : _dof_to_previous_gap[dof_object] = gap;
295 : }
296 78 : }
297 :
298 : void
299 78 : PenaltyWeightedGapUserObject::updateAugmentedLagrangianMultipliers()
300 : {
301 : using std::abs;
302 :
303 3726 : for (const auto & [dof_object, wgap] : _dof_to_weighted_gap)
304 : {
305 : auto & penalty = _dof_to_local_penalty[dof_object];
306 3648 : if (penalty == 0.0)
307 787 : penalty = _penalty;
308 :
309 3648 : 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 3648 : (_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 3648 : if (lagrange_multiplier + gap * penalty / possible_normalization <= 0)
319 282 : lagrange_multiplier += gap * penalty / possible_normalization;
320 : else
321 3366 : lagrange_multiplier = 0.0;
322 :
323 3648 : const auto previous_gap = _dof_to_previous_gap[dof_object];
324 3648 : Real eval_tn = 0;
325 :
326 3648 : 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 3648 : else if (_adaptivity_normal == AdaptivityNormalPenalty::BUSSETTA)
332 3648 : bussettaAdaptivePenalty(previous_gap, gap, penalty, eval_tn);
333 :
334 : // Heuristics to bound the penalty factor
335 3648 : if (penalty < _penalty)
336 0 : penalty = _penalty;
337 3648 : else if (penalty > _max_penalty_multiplier * _penalty)
338 0 : penalty = _max_penalty_multiplier * _penalty;
339 : }
340 78 : }
341 :
342 : void
343 3648 : 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 3648 : if (previous_gap > 0.0)
350 : {
351 2698 : penalty = _penalty;
352 2698 : eval_tn = 0.0;
353 : }
354 : else
355 : {
356 950 : if (previous_gap * gap < 0)
357 22 : eval_tn = 0.0;
358 : else
359 928 : eval_tn = penalty * previous_gap;
360 :
361 950 : adaptiveNormalPenalty(previous_gap, gap, penalty);
362 : }
363 3648 : }
364 :
365 : void
366 950 : PenaltyWeightedGapUserObject::adaptiveNormalPenalty(const Real previous_gap,
367 : const Real gap,
368 : Real & penalty)
369 : {
370 : using std::abs, std::max, std::sqrt;
371 :
372 950 : const bool condition_one = abs(abs(previous_gap / gap) - 1.0) < 0.01;
373 950 : const bool condition_two = abs(gap) > 1.01 * abs(previous_gap);
374 :
375 950 : if (previous_gap * gap < 0)
376 : {
377 22 : 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 22 : penalty = abs(penalty * previous_gap / 10.0 / gap);
382 : }
383 928 : else if (abs(gap) > _penetration_tolerance)
384 : {
385 787 : if (abs(gap - previous_gap) >
386 1635 : max(gap / 10.0, max(previous_gap / 10.0, 5.0 * _penetration_tolerance)))
387 787 : 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 950 : }
|
