www.mooseframework.org
ElementJacobianDamper.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "ElementJacobianDamper.h"
16 #include "FEProblem.h"
17 #include "MooseMesh.h"
18 #include "DisplacedProblem.h"
19 #include "Assembly.h"
20 
21 #include "libmesh/quadrature.h" // _qrule
22 
23 template <>
24 InputParameters
26 {
27  InputParameters params = validParams<GeneralDamper>();
28  params.addClassDescription("Damper that limits the change in element Jacobians");
29  params.addParam<std::vector<NonlinearVariableName>>("displacements",
30  "The nonlinear displacement variables");
31  params.addParam<Real>(
32  "max_increment",
33  0.1,
34  "The maximum permissible relative increment in the Jacobian per Newton iteration");
35  params.set<bool>("use_displaced_mesh") = true;
36  return params;
37 }
38 
39 ElementJacobianDamper::ElementJacobianDamper(const InputParameters & parameters)
40  : GeneralDamper(parameters),
41  _tid(parameters.get<THREAD_ID>("_tid")),
42  _assembly(_subproblem.assembly(_tid)),
43  _qrule(_assembly.qRule()),
44  _JxW(_assembly.JxW()),
45  _fe_problem(*parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
46  _displaced_problem(_fe_problem.getDisplacedProblem()),
47  _max_jacobian_diff(parameters.get<Real>("max_increment"))
48 {
49  if (_displaced_problem == NULL)
50  mooseError("ElementJacobianDamper: Must use displaced problem");
51 
52  _mesh = &_displaced_problem->mesh();
53 
54  const std::vector<NonlinearVariableName> & nl_vnames(
55  getParam<std::vector<NonlinearVariableName>>("displacements"));
56  _ndisp = nl_vnames.size();
57 
58  for (unsigned int i = 0; i < _ndisp; ++i)
59  {
60  _disp_var.push_back(&_sys.getVariable(_tid, nl_vnames[i]));
61  _disp_incr.push_back(_disp_var.back()->increment());
62  }
63 }
64 
65 void
67 {
68 }
69 
70 Real
71 ElementJacobianDamper::computeDamping(const NumericVector<Number> & /* solution */,
72  const NumericVector<Number> & update)
73 {
74  // Maximum difference in the Jacobian for this Newton iteration
75  Real max_difference = 0.0;
76  MooseArray<Real> JxW_displaced;
77 
78  // Vector for storing the original node coordinates
79  std::vector<Point> point_copies;
80 
81  // Loop over elements in the mesh
82  const MeshBase::element_iterator end = _mesh->getMesh().active_local_elements_end();
83  for (auto el = _mesh->getMesh().active_local_elements_begin(); el != end; ++el)
84  {
85  Elem * current_elem = *el;
86 
87  point_copies.clear();
88  point_copies.reserve(current_elem->n_nodes());
89 
90  // Displace nodes with current Newton increment
91  for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
92  {
93  Node & displaced_node = current_elem->node_ref(i);
94 
95  point_copies.push_back(displaced_node);
96 
97  for (unsigned int j = 0; j < _ndisp; ++j)
98  {
99  dof_id_type disp_dof_num =
100  displaced_node.dof_number(_sys.number(), _disp_var[j]->number(), 0);
101  displaced_node(j) += update(disp_dof_num);
102  }
103  }
104 
105  // Reinit element to compute Jacobian of displaced element
106  _assembly.reinit(current_elem);
107  JxW_displaced = _JxW;
108 
109  // Un-displace nodes
110  for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
111  {
112  Node & displaced_node = current_elem->node_ref(i);
113 
114  for (unsigned int j = 0; j < _ndisp; ++j)
115  displaced_node(j) = point_copies[i](j);
116  }
117 
118  // Reinit element to compute Jacobian before displacement
119  _assembly.reinit(current_elem);
120 
121  for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
122  {
123  Real diff = std::abs(JxW_displaced[qp] - _JxW[qp]) / _JxW[qp];
124  if (diff > max_difference)
125  max_difference = diff;
126  }
127 
128  JxW_displaced.release();
129  }
130 
131  _communicator.max(max_difference);
132 
133  if (max_difference > _max_jacobian_diff)
134  return _max_jacobian_diff / max_difference;
135 
136  return 1.0;
137 }
QBase *& _qrule
Quadrature rule.
InputParameters validParams< ElementJacobianDamper >()
unsigned int _ndisp
The number of displacement variables.
MooseSharedPointer< DisplacedProblem > _displaced_problem
The displaced problem.
std::vector< VariableValue > _disp_incr
The current Newton increment in the displacement variables.
MooseMesh * _mesh
The displaced mesh.
const Real _max_jacobian_diff
Maximum allowed relative increment in Jacobian.
const MooseArray< Real > & _JxW
Transformed Jacobian weights.
ElementJacobianDamper(const InputParameters &parameters)
virtual Real computeDamping(const NumericVector< Number > &, const NumericVector< Number > &update) override
Computes this Damper&#39;s damping.
std::vector< MooseVariable * > _disp_var
The displacement variables.
THREAD_ID _tid
Thread ID.
virtual void initialSetup() override