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