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 "ComputeDynamicWeightedGapLMMechanicalContact.h"
11 : #include "MortarContactUtils.h"
12 : #include "DisplacedProblem.h"
13 : #include "Assembly.h"
14 :
15 : #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
16 : #include "metaphysicl/parallel_dualnumber.h"
17 : #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
18 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
19 : #include "timpi/parallel_sync.h"
20 :
21 : registerMooseObject("ContactApp", ComputeDynamicWeightedGapLMMechanicalContact);
22 :
23 : namespace
24 : {
25 : const InputParameters &
26 133 : assignVarsInParamsDyn(const InputParameters & params_in)
27 : {
28 : InputParameters & ret = const_cast<InputParameters &>(params_in);
29 133 : const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
30 133 : if (disp_x_name.size() != 1)
31 0 : mooseError("We require that the disp_x parameter have exactly one coupled name");
32 :
33 : // We do this so we don't get any variable errors during MortarConstraint(Base) construction
34 266 : ret.set<VariableName>("secondary_variable") = disp_x_name[0];
35 133 : ret.set<VariableName>("primary_variable") = disp_x_name[0];
36 :
37 133 : return ret;
38 : }
39 : }
40 : InputParameters
41 266 : ComputeDynamicWeightedGapLMMechanicalContact::validParams()
42 : {
43 266 : InputParameters params = ADMortarConstraint::validParams();
44 266 : params.addClassDescription(
45 : "Computes the normal contact mortar constraints for dynamic simulations");
46 798 : params.addRangeCheckedParam<Real>("capture_tolerance",
47 532 : 1.0e-5,
48 : "capture_tolerance>=0",
49 : "Parameter describing a gap threshold for the application of "
50 : "the persistency constraint in dynamic simulations.");
51 532 : params.addCoupledVar("wear_depth",
52 : "The name of the mortar auxiliary variable that is used to modify the "
53 : "weighted gap definition");
54 532 : params.addRequiredRangeCheckedParam<Real>(
55 : "newmark_beta", "newmark_beta > 0", "Beta parameter for the Newmark time integrator");
56 532 : params.addRequiredRangeCheckedParam<Real>(
57 : "newmark_gamma", "newmark_gamma >= 0.0", "Gamma parameter for the Newmark time integrator");
58 266 : params.suppressParameter<VariableName>("secondary_variable");
59 266 : params.suppressParameter<VariableName>("primary_variable");
60 532 : params.addRequiredCoupledVar("disp_x", "The x displacement variable");
61 532 : params.addRequiredCoupledVar("disp_y", "The y displacement variable");
62 532 : params.addCoupledVar("disp_z", "The z displacement variable");
63 532 : params.addParam<Real>(
64 532 : "c", 1e6, "Parameter for balancing the size of the gap and contact pressure");
65 532 : params.addParam<bool>(
66 : "normalize_c",
67 532 : false,
68 : "Whether to normalize c by weighting function norm. When unnormalized "
69 : "the value of c effectively depends on element size since in the constraint we compare nodal "
70 : "Lagrange Multiplier values to integrated gap values (LM nodal value is independent of "
71 : "element size, where integrated values are dependent on element size).");
72 266 : params.set<bool>("use_displaced_mesh") = true;
73 266 : params.set<bool>("interpolate_normals") = false;
74 266 : return params;
75 0 : }
76 :
77 133 : ComputeDynamicWeightedGapLMMechanicalContact::ComputeDynamicWeightedGapLMMechanicalContact(
78 133 : const InputParameters & parameters)
79 : : ADMortarConstraint(assignVarsInParamsDyn(parameters)),
80 133 : _secondary_disp_x(adCoupledValue("disp_x")),
81 133 : _primary_disp_x(adCoupledNeighborValue("disp_x")),
82 133 : _secondary_disp_y(adCoupledValue("disp_y")),
83 133 : _primary_disp_y(adCoupledNeighborValue("disp_y")),
84 133 : _has_disp_z(isCoupled("disp_z")),
85 133 : _secondary_disp_z(_has_disp_z ? &adCoupledValue("disp_z") : nullptr),
86 133 : _primary_disp_z(_has_disp_z ? &adCoupledNeighborValue("disp_z") : nullptr),
87 266 : _c(getParam<Real>("c")),
88 266 : _normalize_c(getParam<bool>("normalize_c")),
89 133 : _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
90 133 : _disp_x_var(getVar("disp_x", 0)),
91 133 : _disp_y_var(getVar("disp_y", 0)),
92 165 : _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
93 133 : _capture_tolerance(getParam<Real>("capture_tolerance")),
94 133 : _secondary_x_dot(adCoupledDot("disp_x")),
95 133 : _primary_x_dot(adCoupledNeighborValueDot("disp_x")),
96 133 : _secondary_y_dot(adCoupledDot("disp_y")),
97 133 : _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
98 133 : _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
99 133 : _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
100 266 : _has_wear(isParamValid("wear_depth")),
101 170 : _wear_depth(_has_wear ? coupledValueLower("wear_depth") : _zero),
102 266 : _newmark_beta(getParam<Real>("newmark_beta")),
103 532 : _newmark_gamma(getParam<Real>("newmark_gamma"))
104 : {
105 133 : if (!useDual())
106 0 : mooseError("Dynamic mortar contact constraints requires the use of Lagrange multipliers dual "
107 : "interpolation");
108 133 : }
109 :
110 : void
111 749732 : ComputeDynamicWeightedGapLMMechanicalContact::computeQpProperties()
112 : {
113 : // Trim interior node variable derivatives
114 : const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
115 749732 : *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
116 : const auto & secondary_ip_lowerd_map =
117 749732 : amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
118 :
119 749732 : std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
120 : std::array<ADReal, 3> primary_disp{
121 749732 : {_primary_disp_x[_qp], _primary_disp_y[_qp], _has_disp_z ? (*_primary_disp_z)[_qp] : 0}};
122 749732 : std::array<ADReal, 3> secondary_disp{{_secondary_disp_x[_qp],
123 749732 : _secondary_disp_y[_qp],
124 749732 : _has_disp_z ? (*_secondary_disp_z)[_qp] : 0}};
125 :
126 749732 : trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp, false);
127 749732 : trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp, true);
128 :
129 : const ADReal & prim_x = primary_disp[0];
130 : const ADReal & prim_y = primary_disp[1];
131 : const ADReal * prim_z = nullptr;
132 749732 : if (_has_disp_z)
133 : prim_z = &primary_disp[2];
134 :
135 : const ADReal & sec_x = secondary_disp[0];
136 : const ADReal & sec_y = secondary_disp[1];
137 : const ADReal * sec_z = nullptr;
138 749732 : if (_has_disp_z)
139 : sec_z = &secondary_disp[2];
140 :
141 : std::array<ADReal, 3> primary_disp_dot{
142 749732 : {_primary_x_dot[_qp], _primary_y_dot[_qp], _has_disp_z ? (*_primary_z_dot)[_qp] : 0}};
143 : std::array<ADReal, 3> secondary_disp_dot{
144 749732 : {_secondary_x_dot[_qp], _secondary_y_dot[_qp], _has_disp_z ? (*_secondary_z_dot)[_qp] : 0}};
145 :
146 749732 : trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp_dot, false);
147 749732 : trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp_dot, true);
148 :
149 : const ADReal & prim_x_dot = primary_disp_dot[0];
150 : const ADReal & prim_y_dot = primary_disp_dot[1];
151 : const ADReal * prim_z_dot = nullptr;
152 749732 : if (_has_disp_z)
153 : prim_z_dot = &primary_disp_dot[2];
154 :
155 : const ADReal & sec_x_dot = secondary_disp_dot[0];
156 : const ADReal & sec_y_dot = secondary_disp_dot[1];
157 : const ADReal * sec_z_dot = nullptr;
158 749732 : if (_has_disp_z)
159 : sec_z_dot = &secondary_disp_dot[2];
160 :
161 : // Compute dynamic constraint-related quantities
162 749732 : ADRealVectorValue gap_vec = _phys_points_primary[_qp] - _phys_points_secondary[_qp];
163 749732 : gap_vec(0).derivatives() = prim_x.derivatives() - sec_x.derivatives();
164 749732 : gap_vec(1).derivatives() = prim_y.derivatives() - sec_y.derivatives();
165 :
166 749732 : _relative_velocity = ADRealVectorValue(prim_x_dot - sec_x_dot, prim_y_dot - sec_y_dot, 0.0);
167 :
168 749732 : if (_has_disp_z)
169 : {
170 537632 : gap_vec(2).derivatives() = prim_z->derivatives() - sec_z->derivatives();
171 537632 : _relative_velocity(2) = *prim_z_dot - *sec_z_dot;
172 : }
173 :
174 749732 : _qp_gap_nodal = gap_vec * (_JxW_msm[_qp] * _coord[_qp]);
175 749732 : _qp_velocity = _relative_velocity * (_JxW_msm[_qp] * _coord[_qp]);
176 :
177 : // Current part of the gap velocity Newmark-beta time discretization
178 : _qp_gap_nodal_dynamics =
179 749732 : (_newmark_gamma / _newmark_beta * gap_vec / _dt) * (_JxW_msm[_qp] * _coord[_qp]);
180 :
181 : // To do normalization of constraint coefficient (c_n)
182 749732 : _qp_factor = _JxW_msm[_qp] * _coord[_qp];
183 749732 : }
184 :
185 : ADReal
186 0 : ComputeDynamicWeightedGapLMMechanicalContact::computeQpResidual(Moose::MortarType)
187 : {
188 0 : mooseError(
189 : "We should never call computeQpResidual for ComputeDynamicWeightedGapLMMechanicalContact");
190 : }
191 :
192 : void
193 2574728 : ComputeDynamicWeightedGapLMMechanicalContact::computeQpIProperties()
194 : {
195 : mooseAssert(_normals.size() == _lower_secondary_elem->n_nodes(),
196 : "Making sure that _normals is the expected size");
197 :
198 : // Get the _dof_to_weighted_gap map
199 2574728 : const DofObject * dof = _var->isNodal()
200 2574728 : ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
201 0 : : static_cast<const DofObject *>(_lower_secondary_elem);
202 :
203 : // Regular normal contact constraint: Use before contact is established for contact detection
204 5149456 : _dof_to_weighted_gap[dof].first += _test[_i][_qp] * (_qp_gap_nodal * _normals[_i]);
205 :
206 : // Integrated part of the "persistency" constraint
207 5149456 : _dof_to_weighted_gap_dynamics[dof] += _test[_i][_qp] * _qp_gap_nodal_dynamics * _normals[_i];
208 5149456 : _dof_to_velocity[dof] += _test[_i][_qp] * _qp_velocity * _normals[_i];
209 :
210 2574728 : _dof_to_nodal_wear_depth[dof] += _test[_i][_qp] * _wear_depth[_qp] * _JxW_msm[_qp] * _coord[_qp];
211 :
212 2574728 : if (_normalize_c)
213 72400 : _dof_to_weighted_gap[dof].second += _test[_i][_qp] * _qp_factor;
214 2574728 : }
215 :
216 : void
217 568 : ComputeDynamicWeightedGapLMMechanicalContact::timestepSetup()
218 : {
219 : // These dof maps are not recoverable as they are maps of pointers, and recovering old pointers
220 : // would be wrong. We would need to create a custom dataStore() and dataLoad()
221 568 : if (_app.isRecovering())
222 0 : mooseError("This object does not support recovering");
223 :
224 : _dof_to_old_weighted_gap.clear();
225 : _dof_to_old_velocity.clear();
226 : _dof_to_nodal_old_wear_depth.clear();
227 :
228 3583 : for (auto & map_pr : _dof_to_weighted_gap)
229 3015 : _dof_to_old_weighted_gap.emplace(map_pr.first, std::move(map_pr.second.first));
230 :
231 3583 : for (auto & map_pr : _dof_to_velocity)
232 : _dof_to_old_velocity.emplace(map_pr);
233 :
234 3583 : for (auto & map_pr : _dof_to_nodal_wear_depth)
235 : _dof_to_nodal_old_wear_depth.emplace(map_pr);
236 568 : }
237 :
238 : void
239 10020 : ComputeDynamicWeightedGapLMMechanicalContact::residualSetup()
240 : {
241 : _dof_to_weighted_gap.clear();
242 : _dof_to_weighted_gap_dynamics.clear();
243 : _dof_to_velocity.clear();
244 :
245 : // Wear
246 : _dof_to_nodal_wear_depth.clear();
247 10020 : }
248 :
249 : void
250 608 : ComputeDynamicWeightedGapLMMechanicalContact::post()
251 : {
252 608 : Moose::Mortar::Contact::communicateGaps(
253 608 : _dof_to_weighted_gap, _mesh, _nodal, _normalize_c, _communicator, false);
254 :
255 608 : if (_has_wear)
256 0 : communicateWear();
257 :
258 : // There is a need for the dynamic constraint to uncouple the computation of the weighted gap from
259 : // the computation of the constraint itself since we are switching from gap constraint to
260 : // persistency constraint.
261 4028 : for (const auto & pr : _dof_to_weighted_gap)
262 : {
263 3420 : if (pr.first->processor_id() != this->processor_id())
264 : continue;
265 :
266 : //
267 6840 : _dof_to_weighted_gap[pr.first].first += _dof_to_nodal_wear_depth[pr.first];
268 : _dof_to_weighted_gap_dynamics[pr.first] +=
269 6840 : _newmark_gamma / _newmark_beta * _dof_to_nodal_wear_depth[pr.first] / _dt;
270 : //
271 :
272 : const auto is_dof_on_map = _dof_to_old_weighted_gap.find(pr.first);
273 :
274 : // If is_dof_on_map isn't on map, it means it's an initial step
275 4095 : if (is_dof_on_map == _dof_to_old_weighted_gap.end() ||
276 : _dof_to_old_weighted_gap[pr.first] > _capture_tolerance)
277 2745 : _weighted_gap_ptr = &pr.second.first;
278 : else
279 : {
280 1350 : ADReal term = -_newmark_gamma / _newmark_beta / _dt * _dof_to_old_weighted_gap[pr.first];
281 675 : term += _dof_to_old_velocity[pr.first];
282 675 : _dof_to_weighted_gap_dynamics[pr.first] += term;
283 675 : _weighted_gap_ptr = &_dof_to_weighted_gap_dynamics[pr.first];
284 : }
285 :
286 3420 : _normalization_ptr = &pr.second.second;
287 :
288 3420 : ComputeDynamicWeightedGapLMMechanicalContact::enforceConstraintOnDof(pr.first);
289 : }
290 608 : }
291 :
292 : void
293 9412 : ComputeDynamicWeightedGapLMMechanicalContact::incorrectEdgeDroppingPost(
294 : const std::unordered_set<const Node *> & inactive_lm_nodes)
295 : {
296 9412 : Moose::Mortar::Contact::communicateGaps(
297 9412 : _dof_to_weighted_gap, _mesh, _nodal, _normalize_c, _communicator, false);
298 :
299 9412 : if (_has_wear)
300 1448 : communicateWear();
301 :
302 77017 : for (const auto & pr : _dof_to_weighted_gap)
303 : {
304 67605 : if ((inactive_lm_nodes.find(static_cast<const Node *>(pr.first)) != inactive_lm_nodes.end()) ||
305 : (pr.first->processor_id() != this->processor_id()))
306 : continue;
307 :
308 : //
309 132946 : _dof_to_weighted_gap[pr.first].first += _dof_to_nodal_wear_depth[pr.first];
310 : _dof_to_weighted_gap_dynamics[pr.first] +=
311 132946 : _newmark_gamma / _newmark_beta * _dof_to_nodal_wear_depth[pr.first] / _dt;
312 : const auto is_dof_on_map = _dof_to_old_weighted_gap.find(pr.first);
313 :
314 : // If is_dof_on_map isn't on map, it means it's an initial step
315 116253 : if (is_dof_on_map == _dof_to_old_weighted_gap.end() ||
316 : _dof_to_old_weighted_gap[pr.first] > _capture_tolerance)
317 : {
318 : // If this is the first step or the previous step gap is not identified as in contact, apply
319 : // regular conditions
320 50633 : _weighted_gap_ptr = &pr.second.first;
321 : }
322 : else
323 : {
324 15840 : ADReal term = _dof_to_weighted_gap[pr.first].first * _newmark_gamma / (_newmark_beta * _dt);
325 31680 : term -= _dof_to_old_weighted_gap[pr.first] * _newmark_gamma / (_newmark_beta * _dt);
326 15840 : term -= _dof_to_old_velocity[pr.first];
327 15840 : _dof_to_weighted_gap_dynamics[pr.first] = term;
328 : // Enable the application of persistency condition
329 15840 : _weighted_gap_ptr = &_dof_to_weighted_gap_dynamics[pr.first];
330 : }
331 :
332 66473 : _normalization_ptr = &pr.second.second;
333 :
334 66473 : ComputeDynamicWeightedGapLMMechanicalContact::enforceConstraintOnDof(pr.first);
335 : }
336 9412 : }
337 :
338 : void
339 1448 : ComputeDynamicWeightedGapLMMechanicalContact::communicateWear()
340 : {
341 : // We may have wear depth information that should go to other processes that own the dofs
342 : using Datum = std::pair<dof_id_type, ADReal>;
343 : std::unordered_map<processor_id_type, std::vector<Datum>> push_data;
344 :
345 11403 : for (auto & pr : _dof_to_nodal_wear_depth)
346 : {
347 9955 : const auto * const dof_object = pr.first;
348 9955 : const auto proc_id = dof_object->processor_id();
349 9955 : if (proc_id == this->processor_id())
350 9955 : continue;
351 :
352 0 : push_data[proc_id].push_back(std::make_pair(dof_object->id(), std::move(pr.second)));
353 : }
354 :
355 1448 : const auto & lm_mesh = _mesh.getMesh();
356 :
357 0 : auto action_functor = [this, &lm_mesh](const processor_id_type libmesh_dbg_var(pid),
358 0 : const std::vector<Datum> & sent_data)
359 : {
360 : mooseAssert(pid != this->processor_id(), "We do not send messages to ourself here");
361 0 : for (auto & pr : sent_data)
362 : {
363 0 : const auto dof_id = pr.first;
364 : const auto * const dof_object =
365 0 : _nodal ? static_cast<const DofObject *>(lm_mesh.node_ptr(dof_id))
366 0 : : static_cast<const DofObject *>(lm_mesh.elem_ptr(dof_id));
367 : mooseAssert(dof_object, "This should be non-null");
368 0 : _dof_to_nodal_wear_depth[dof_object] += std::move(pr.second);
369 : }
370 1448 : };
371 :
372 1448 : TIMPI::push_parallel_vector_data(_communicator, push_data, action_functor);
373 1448 : }
374 :
375 : void
376 69893 : ComputeDynamicWeightedGapLMMechanicalContact::enforceConstraintOnDof(const DofObject * const dof)
377 : {
378 69893 : const auto & weighted_gap = *_weighted_gap_ptr;
379 69893 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
380 :
381 69893 : const auto dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
382 69893 : ADReal lm_value = (*_sys.currentSolution())(dof_index);
383 : Moose::derivInsert(lm_value.derivatives(), dof_index, 1.);
384 :
385 69893 : const ADReal dof_residual = std::min(lm_value, weighted_gap * c);
386 :
387 69893 : addResidualsAndJacobian(_assembly,
388 69893 : std::array<ADReal, 1>{{dof_residual}},
389 139786 : std::array<dof_id_type, 1>{{dof_index}},
390 69893 : _var->scalingFactor());
391 69893 : }
392 :
393 : void
394 684574 : ComputeDynamicWeightedGapLMMechanicalContact::computeResidual(const Moose::MortarType mortar_type)
395 : {
396 684574 : if (mortar_type != Moose::MortarType::Lower)
397 : return;
398 :
399 : mooseAssert(_var, "LM variable is null");
400 :
401 990190 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
402 : {
403 749732 : computeQpProperties();
404 3324460 : for (_i = 0; _i < _test.size(); ++_i)
405 2574728 : computeQpIProperties();
406 : }
407 : }
408 :
409 : void
410 2780 : ComputeDynamicWeightedGapLMMechanicalContact::jacobianSetup()
411 : {
412 2780 : residualSetup();
413 2780 : }
414 :
415 : void
416 204323 : ComputeDynamicWeightedGapLMMechanicalContact::computeJacobian(const Moose::MortarType mortar_type)
417 : {
418 : // During "computeResidual" and "computeJacobian" we are actually just computing properties on the
419 : // mortar segment element mesh. We are *not* actually assembling into the residual/Jacobian. For
420 : // the zero-penetration constraint, the property of interest is the map from node to weighted gap.
421 : // Computation of the properties proceeds identically for residual and Jacobian evaluation hence
422 : // why we simply call computeResidual here. We will assemble into the residual/Jacobian later from
423 : // the post() method
424 204323 : computeResidual(mortar_type);
425 204323 : }
|
