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 : // compute new normal pressure for each node
136 146402 : for (const auto & [dof_object, wgap] : _dof_to_weighted_gap)
137 : {
138 116282 : 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 116282 : auto gap = _use_physical_gap ? adPhysicalGap(wgap) / wgap.second : wgap.first;
142 :
143 : const auto lagrange_multiplier =
144 176154 : _augmented_lagrange_problem ? _dof_to_lagrange_multiplier[dof_object] : 0.0;
145 :
146 : // keep the negative normal pressure (compressive per convention here)
147 116282 : 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 232564 : _dof_to_normal_pressure[dof_object] = -normal_pressure;
151 : }
152 30120 : }
153 :
154 : void
155 13424 : PenaltyWeightedGapUserObject::finalize()
156 : {
157 13424 : WeightedGapUserObject::finalize();
158 13424 : selfFinalize();
159 13424 : }
160 :
161 : Real
162 157098 : PenaltyWeightedGapUserObject::getNormalContactPressure(const Node * const node) const
163 : {
164 157098 : const auto it = _dof_to_normal_pressure.find(_subproblem.mesh().nodePtr(node->id()));
165 :
166 157098 : if (it != _dof_to_normal_pressure.end())
167 4496 : return MetaPhysicL::raw_value(it->second);
168 : else
169 : return 0.0;
170 : }
171 :
172 : Real
173 972 : PenaltyWeightedGapUserObject::getNormalLagrangeMultiplier(const Node * const node) const
174 : {
175 972 : const auto it = _dof_to_lagrange_multiplier.find(_subproblem.mesh().nodePtr(node->id()));
176 :
177 972 : if (it != _dof_to_lagrange_multiplier.end())
178 852 : return -MetaPhysicL::raw_value(it->second);
179 : else
180 : return 0.0;
181 : }
182 :
183 : void
184 206880 : PenaltyWeightedGapUserObject::reinit()
185 : {
186 206880 : _contact_pressure.resize(_qrule_msm->n_points());
187 778336 : for (const auto qp : make_range(_qrule_msm->n_points()))
188 571456 : _contact_pressure[qp] = 0.0;
189 :
190 778336 : for (const auto i : make_range(_test->size()))
191 : {
192 571456 : const Node * const node = _lower_secondary_elem->node_ptr(i);
193 :
194 2345152 : for (const auto qp : make_range(_qrule_msm->n_points()))
195 3547392 : _contact_pressure[qp] += (*_test)[i][qp] * _dof_to_normal_pressure[node];
196 : }
197 206880 : }
198 :
199 : void
200 2062 : PenaltyWeightedGapUserObject::selfTimestepSetup()
201 : {
202 : // Let's not clear the LMs for improved performance in
203 : // nonlinear problems
204 :
205 : // reset penalty
206 2133 : for (auto & dof_lp : _dof_to_local_penalty)
207 71 : dof_lp.second = _penalty;
208 :
209 : // clear previous gap
210 2133 : for (auto & dof_pg : _dof_to_previous_gap)
211 71 : dof_pg.second = 0.0;
212 :
213 : // save old timestep
214 2062 : _dt_old = _dt;
215 2062 : }
216 :
217 : void
218 1067 : PenaltyWeightedGapUserObject::timestepSetup()
219 : {
220 1067 : selfTimestepSetup();
221 1067 : }
222 :
223 : bool
224 96 : PenaltyWeightedGapUserObject::isAugmentedLagrangianConverged()
225 : {
226 96 : Real max_positive_gap = 0.0;
227 96 : Real min_negative_gap = 0.0;
228 :
229 : // Get maximum gap to ascertain whether we are converged.
230 4602 : for (auto & [dof_object, wgap] : _dof_to_weighted_gap)
231 : {
232 4506 : const auto gap = physicalGap(wgap);
233 : {
234 : // Check condition for nodes that are active
235 8768 : if (gap < 0 || _dof_to_lagrange_multiplier[dof_object] < 0.0)
236 : {
237 348 : max_positive_gap = std::max(max_positive_gap, gap);
238 348 : 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 96 : if (this->_communicator.rank() == 0)
246 63 : recv.resize(this->_communicator.size());
247 96 : this->_communicator.gather(
248 96 : 0, -min_negative_gap > max_positive_gap ? min_negative_gap : max_positive_gap, recv);
249 :
250 96 : if (this->_communicator.rank() == 0)
251 : {
252 63 : min_negative_gap = *std::min_element(recv.begin(), recv.end());
253 63 : max_positive_gap = *std::max_element(recv.begin(), recv.end());
254 :
255 : // report the gap value with the largest magnitude
256 63 : const Real reported_gap =
257 63 : -min_negative_gap > max_positive_gap ? min_negative_gap : max_positive_gap;
258 63 : if (std::abs(reported_gap) > _penetration_tolerance)
259 : {
260 11 : mooseInfoRepeated("Penetration tolerance fail max_gap = ",
261 : reported_gap,
262 : " (gap_tol=",
263 11 : _penetration_tolerance,
264 : "). Iteration number is: ",
265 : _lagrangian_iteration_number,
266 : ".");
267 11 : return false;
268 : }
269 : else
270 52 : mooseInfoRepeated("Penetration tolerance success max_gap = ",
271 : reported_gap,
272 : " (gap_tol=",
273 52 : _penetration_tolerance,
274 : "). Iteration number is: ",
275 : _lagrangian_iteration_number,
276 : ".");
277 : }
278 :
279 : return true;
280 96 : }
281 :
282 : void
283 78 : PenaltyWeightedGapUserObject::augmentedLagrangianSetup()
284 : {
285 : // Loop over all nodes for which a gap has been computed
286 3726 : for (auto & [dof_object, wgap] : _dof_to_weighted_gap)
287 : {
288 : const Real gap = physicalGap(wgap);
289 : // Store previous augmented lagrange iteration gap
290 3648 : _dof_to_previous_gap[dof_object] = gap;
291 : }
292 78 : }
293 :
294 : void
295 78 : PenaltyWeightedGapUserObject::updateAugmentedLagrangianMultipliers()
296 : {
297 3726 : for (const auto & [dof_object, wgap] : _dof_to_weighted_gap)
298 : {
299 : auto & penalty = _dof_to_local_penalty[dof_object];
300 3648 : if (penalty == 0.0)
301 787 : penalty = _penalty;
302 :
303 3648 : 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 3648 : (_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 3648 : if (lagrange_multiplier + gap * penalty / possible_normalization <= 0)
313 282 : lagrange_multiplier += gap * penalty / possible_normalization;
314 : else
315 3366 : lagrange_multiplier = 0.0;
316 :
317 3648 : const auto previous_gap = _dof_to_previous_gap[dof_object];
318 3648 : Real eval_tn = 0;
319 :
320 3648 : 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 3648 : else if (_adaptivity_normal == AdaptivityNormalPenalty::BUSSETTA)
326 3648 : bussettaAdaptivePenalty(previous_gap, gap, penalty, eval_tn);
327 :
328 : // Heuristics to bound the penalty factor
329 3648 : if (penalty < _penalty)
330 0 : penalty = _penalty;
331 3648 : else if (penalty > _max_penalty_multiplier * _penalty)
332 0 : penalty = _max_penalty_multiplier * _penalty;
333 : }
334 78 : }
335 :
336 : void
337 3648 : 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 3648 : if (previous_gap > 0.0)
344 : {
345 2698 : penalty = _penalty;
346 2698 : eval_tn = 0.0;
347 : }
348 : else
349 : {
350 950 : if (previous_gap * gap < 0)
351 22 : eval_tn = 0.0;
352 : else
353 928 : eval_tn = penalty * previous_gap;
354 :
355 950 : adaptiveNormalPenalty(previous_gap, gap, penalty);
356 : }
357 3648 : }
358 :
359 : void
360 950 : PenaltyWeightedGapUserObject::adaptiveNormalPenalty(const Real previous_gap,
361 : const Real gap,
362 : Real & penalty)
363 : {
364 950 : const bool condition_one = std::abs(std::abs(previous_gap / gap) - 1.0) < 0.01;
365 950 : const bool condition_two = std::abs(gap) > 1.01 * std::abs(previous_gap);
366 :
367 950 : if (previous_gap * gap < 0)
368 : {
369 22 : 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 22 : penalty = std::abs(penalty * previous_gap / 10.0 / gap);
374 : }
375 928 : else if (std::abs(gap) > _penetration_tolerance)
376 : {
377 787 : if (std::abs(gap - previous_gap) >
378 1635 : std::max(gap / 10.0, std::max(previous_gap / 10.0, 5.0 * _penetration_tolerance)))
379 787 : 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 950 : }
|
