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 "ElementJacobianDamper.h"
11 : #include "FEProblem.h"
12 : #include "MooseMesh.h"
13 : #include "DisplacedProblem.h"
14 : #include "Assembly.h"
15 :
16 : #include "libmesh/quadrature.h" // _qrule
17 :
18 : // C++
19 : #include <algorithm>
20 : #include <cstring> // for "Jacobian" exception test
21 : #include <limits>
22 :
23 : registerMooseObject("SolidMechanicsApp", ElementJacobianDamper);
24 :
25 : namespace
26 : {
27 : void
28 12829 : restoreNodes(Elem & elem, const std::vector<Point> & point_copies)
29 : {
30 : mooseAssert(elem.n_nodes() == point_copies.size(), "Node restore cache is the wrong size");
31 :
32 115461 : for (const auto i : make_range(elem.n_nodes()))
33 102632 : elem.node_ref(i) = point_copies[i];
34 12829 : }
35 : }
36 :
37 : InputParameters
38 22 : ElementJacobianDamper::validParams()
39 : {
40 22 : InputParameters params = GeneralDamper::validParams();
41 22 : params.addClassDescription("Damper that limits the change in element Jacobians");
42 44 : params.addParam<std::vector<VariableName>>(
43 : "displacements", {}, "The nonlinear displacement variables");
44 44 : params.addParam<Real>(
45 : "max_increment",
46 44 : 0.1,
47 : "The maximum permissible relative increment in the Jacobian per Newton iteration");
48 44 : params.addParam<bool>(
49 : "use_backtracking",
50 44 : false,
51 : "If true, iteratively cut back the probed Newton update when a full trial update would "
52 : "create a degenerate displaced element map. The accepted trial must also satisfy "
53 : "max_increment.");
54 66 : params.addRangeCheckedParam<Real>(
55 : "backtrack_factor",
56 44 : 0.5,
57 : "backtrack_factor > 0 & backtrack_factor < 1",
58 : "Multiplicative cutback applied to the trial damping during ElementJacobianDamper "
59 : "backtracking.");
60 44 : params.addParam<unsigned int>(
61 44 : "max_backtrack_steps", 12, "Maximum number of ElementJacobianDamper backtracking attempts.");
62 22 : params.set<bool>("use_displaced_mesh") = true;
63 22 : return params;
64 0 : }
65 :
66 11 : ElementJacobianDamper::ElementJacobianDamper(const InputParameters & parameters)
67 : : GeneralDamper(parameters),
68 11 : _tid(parameters.get<THREAD_ID>("_tid")),
69 11 : _assembly(_subproblem.assembly(_tid, _sys.number())),
70 11 : _qrule(_assembly.qRule()),
71 11 : _JxW(_assembly.JxW()),
72 11 : _fe_problem(*parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
73 11 : _displaced_problem(_fe_problem.getDisplacedProblem()),
74 22 : _max_jacobian_diff(getParam<Real>("max_increment")),
75 22 : _use_backtracking(getParam<bool>("use_backtracking")),
76 22 : _backtrack_factor(getParam<Real>("backtrack_factor")),
77 33 : _max_backtrack_steps(getParam<unsigned int>("max_backtrack_steps"))
78 : {
79 11 : if (_displaced_problem == NULL)
80 0 : mooseError("ElementJacobianDamper: Must use displaced problem");
81 :
82 11 : _mesh = &_displaced_problem->mesh();
83 :
84 22 : const std::vector<VariableName> & nl_vnames(getParam<std::vector<VariableName>>("displacements"));
85 11 : _ndisp = nl_vnames.size();
86 :
87 44 : for (unsigned int i = 0; i < _ndisp; ++i)
88 : {
89 33 : _disp_var.push_back(&_sys.getFieldVariable<Real>(_tid, nl_vnames[i]));
90 33 : _disp_incr.push_back(_disp_var.back()->increment());
91 : }
92 11 : }
93 :
94 : void
95 11 : ElementJacobianDamper::initialSetup()
96 : {
97 11 : }
98 :
99 : Real
100 243 : ElementJacobianDamper::computeDamping(const NumericVector<Number> & /* solution */,
101 : const NumericVector<Number> & update)
102 : {
103 : auto probe_trial =
104 286 : [&](const Real trial_damping, Real & max_difference, std::string & error_message)
105 : {
106 : const bool trial_nondegenerate_local =
107 286 : probeDamping(update, trial_damping, max_difference, error_message);
108 :
109 286 : _communicator.max(max_difference);
110 :
111 286 : const unsigned int invalid_local = trial_nondegenerate_local ? 0 : 1;
112 286 : unsigned int invalid = invalid_local;
113 286 : _communicator.max(invalid);
114 :
115 286 : if (invalid)
116 : {
117 : processor_id_type first_invalid_processor =
118 8 : invalid_local ? _communicator.rank() : std::numeric_limits<processor_id_type>::max();
119 8 : _communicator.min(first_invalid_processor);
120 8 : _communicator.broadcast(error_message, first_invalid_processor);
121 : }
122 :
123 286 : return !invalid;
124 243 : };
125 :
126 243 : if (!_use_backtracking)
127 : {
128 218 : Real max_difference = 0.0;
129 : std::string error_message;
130 :
131 : PARALLEL_TRY
132 : {
133 218 : if (!probe_trial(/*trial_damping=*/1.0, max_difference, error_message))
134 3 : _fe_problem.setException(error_message);
135 : }
136 218 : PARALLEL_CATCH;
137 :
138 215 : if (max_difference > _max_jacobian_diff)
139 159 : return _max_jacobian_diff / max_difference;
140 :
141 : return 1.0;
142 : }
143 :
144 50 : const Real minimum_trial_damping = std::max(_min_damping, std::numeric_limits<Real>::epsilon());
145 : Real damping = 1.0;
146 : std::string error_message;
147 :
148 : // Try 1 step, then at most max_backtrack_steps extra steps
149 68 : for (unsigned int step = 0; step <= _max_backtrack_steps; ++step)
150 : {
151 68 : Real max_difference = 0.0;
152 : bool trial_nondegenerate = false;
153 :
154 68 : PARALLEL_TRY { trial_nondegenerate = probe_trial(damping, max_difference, error_message); }
155 68 : PARALLEL_CATCH;
156 :
157 68 : if (trial_nondegenerate && max_difference <= _max_jacobian_diff)
158 25 : return damping;
159 :
160 43 : if (damping <= minimum_trial_damping || step == _max_backtrack_steps)
161 : {
162 0 : if (!trial_nondegenerate)
163 0 : _fe_problem.setException(error_message.empty()
164 0 : ? "ElementJacobianDamper could not find a nondegenerate "
165 : "trial update."
166 : : error_message);
167 : else
168 0 : _fe_problem.setException("ElementJacobianDamper could not reduce the relative Jacobian "
169 : "increment below max_increment without driving the damping "
170 : "factor below a usable threshold.");
171 :
172 : // The exception flag will abort the nonlinear step after the damper returns, so return the
173 : // last trial damping only to satisfy the interface.
174 0 : return damping;
175 : }
176 :
177 43 : if (!trial_nondegenerate)
178 5 : damping *= _backtrack_factor;
179 : else
180 38 : damping =
181 61 : std::min(damping * _backtrack_factor, damping * _max_jacobian_diff / max_difference);
182 : }
183 :
184 : return damping;
185 : }
186 :
187 : bool
188 286 : ElementJacobianDamper::probeDamping(const NumericVector<Number> & update,
189 : const Real damping,
190 : Real & max_difference,
191 : std::string & error_message)
192 : {
193 : MooseArray<Real> JxW_displaced;
194 :
195 : // Vector for storing the original node coordinates
196 : std::vector<Point> point_copies;
197 :
198 286 : max_difference = 0.0;
199 : error_message.clear();
200 :
201 : try
202 : {
203 : // Loop over elements in the mesh
204 26216 : for (auto & current_elem : _mesh->getMesh().active_local_element_ptr_range())
205 : {
206 12829 : point_copies.clear();
207 12829 : point_copies.reserve(current_elem->n_nodes());
208 :
209 : // Displace nodes with the trial Newton increment
210 115461 : for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
211 : {
212 102632 : Node & displaced_node = current_elem->node_ref(i);
213 :
214 102632 : point_copies.push_back(displaced_node);
215 :
216 410528 : for (unsigned int j = 0; j < _ndisp; ++j)
217 : {
218 : const dof_id_type disp_dof_num =
219 307896 : displaced_node.dof_number(_sys.number(), _disp_var[j]->number(), 0);
220 307896 : displaced_node(j) += damping * update(disp_dof_num);
221 : }
222 : }
223 :
224 : try
225 : {
226 : // Reinit element to compute Jacobian of the trial displaced element
227 12829 : _assembly.reinit(current_elem);
228 12822 : JxW_displaced = _JxW;
229 : }
230 7 : catch (const std::exception & e)
231 : {
232 7 : restoreNodes(*current_elem, point_copies);
233 :
234 : // Degenerate-map failures mean this trial damping is too aggressive
235 7 : if (strstr(e.what(), "Jacobian") || strstr(e.what(), "singular") ||
236 0 : strstr(e.what(), "det != 0"))
237 : {
238 7 : error_message = e.what();
239 : return false;
240 : }
241 :
242 0 : throw;
243 7 : }
244 :
245 12822 : restoreNodes(*current_elem, point_copies);
246 :
247 : // Reinit element to compute Jacobian before displacement
248 : try
249 : {
250 12822 : _assembly.reinit(current_elem);
251 : }
252 0 : catch (const std::exception & e)
253 : {
254 0 : if (strstr(e.what(), "Jacobian") || strstr(e.what(), "singular") ||
255 0 : strstr(e.what(), "det != 0"))
256 : {
257 0 : error_message = e.what();
258 : return false;
259 : }
260 :
261 0 : throw;
262 0 : }
263 :
264 115398 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
265 : {
266 102576 : const Real denominator = std::max(std::abs(_JxW[qp]), libMesh::TOLERANCE);
267 102576 : const Real diff = std::abs(JxW_displaced[qp] - _JxW[qp]) / denominator;
268 102576 : if (diff > max_difference)
269 970 : max_difference = diff;
270 : }
271 :
272 12822 : JxW_displaced.release();
273 286 : }
274 : }
275 0 : catch (const MooseException & e)
276 : {
277 0 : error_message = e.what();
278 0 : _fe_problem.setException(e.what());
279 : return false;
280 0 : }
281 :
282 279 : return true;
283 286 : }
|
