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