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             : registerMooseObject("SolidMechanicsApp", ElementJacobianDamper);
      19             : 
      20             : InputParameters
      21          26 : ElementJacobianDamper::validParams()
      22             : {
      23          26 :   InputParameters params = GeneralDamper::validParams();
      24          26 :   params.addClassDescription("Damper that limits the change in element Jacobians");
      25          52 :   params.addParam<std::vector<VariableName>>(
      26             :       "displacements", {}, "The nonlinear displacement variables");
      27          52 :   params.addParam<Real>(
      28             :       "max_increment",
      29          52 :       0.1,
      30             :       "The maximum permissible relative increment in the Jacobian per Newton iteration");
      31          26 :   params.set<bool>("use_displaced_mesh") = true;
      32          26 :   return params;
      33           0 : }
      34             : 
      35          13 : ElementJacobianDamper::ElementJacobianDamper(const InputParameters & parameters)
      36             :   : GeneralDamper(parameters),
      37          13 :     _tid(parameters.get<THREAD_ID>("_tid")),
      38          13 :     _assembly(_subproblem.assembly(_tid, _sys.number())),
      39          13 :     _qrule(_assembly.qRule()),
      40          13 :     _JxW(_assembly.JxW()),
      41          13 :     _fe_problem(*parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
      42          13 :     _displaced_problem(_fe_problem.getDisplacedProblem()),
      43          26 :     _max_jacobian_diff(parameters.get<Real>("max_increment"))
      44             : {
      45          13 :   if (_displaced_problem == NULL)
      46           0 :     mooseError("ElementJacobianDamper: Must use displaced problem");
      47             : 
      48          13 :   _mesh = &_displaced_problem->mesh();
      49             : 
      50          26 :   const std::vector<VariableName> & nl_vnames(getParam<std::vector<VariableName>>("displacements"));
      51          13 :   _ndisp = nl_vnames.size();
      52             : 
      53          52 :   for (unsigned int i = 0; i < _ndisp; ++i)
      54             :   {
      55          39 :     _disp_var.push_back(&_sys.getFieldVariable<Real>(_tid, nl_vnames[i]));
      56          39 :     _disp_incr.push_back(_disp_var.back()->increment());
      57             :   }
      58          13 : }
      59             : 
      60             : void
      61          13 : ElementJacobianDamper::initialSetup()
      62             : {
      63          13 : }
      64             : 
      65             : Real
      66         284 : ElementJacobianDamper::computeDamping(const NumericVector<Number> & /* solution */,
      67             :                                       const NumericVector<Number> & update)
      68             : {
      69             :   // Maximum difference in the Jacobian for this Newton iteration
      70         284 :   Real max_difference = 0.0;
      71             :   MooseArray<Real> JxW_displaced;
      72             : 
      73             :   // Vector for storing the original node coordinates
      74             :   std::vector<Point> point_copies;
      75             : 
      76             :   // Loop over elements in the mesh
      77         968 :   for (auto & current_elem : _mesh->getMesh().active_local_element_ptr_range())
      78             :   {
      79         200 :     point_copies.clear();
      80         200 :     point_copies.reserve(current_elem->n_nodes());
      81             : 
      82             :     // Displace nodes with current Newton increment
      83        1800 :     for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
      84             :     {
      85        1600 :       Node & displaced_node = current_elem->node_ref(i);
      86             : 
      87        1600 :       point_copies.push_back(displaced_node);
      88             : 
      89        6400 :       for (unsigned int j = 0; j < _ndisp; ++j)
      90             :       {
      91             :         dof_id_type disp_dof_num =
      92        4800 :             displaced_node.dof_number(_sys.number(), _disp_var[j]->number(), 0);
      93        4800 :         displaced_node(j) += update(disp_dof_num);
      94             :       }
      95             :     }
      96             : 
      97             :     // Reinit element to compute Jacobian of displaced element
      98         200 :     _assembly.reinit(current_elem);
      99         200 :     JxW_displaced = _JxW;
     100             : 
     101             :     // Un-displace nodes
     102        1800 :     for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
     103             :     {
     104        1600 :       Node & displaced_node = current_elem->node_ref(i);
     105             : 
     106        6400 :       for (unsigned int j = 0; j < _ndisp; ++j)
     107        4800 :         displaced_node(j) = point_copies[i](j);
     108             :     }
     109             : 
     110             :     // Reinit element to compute Jacobian before displacement
     111         200 :     _assembly.reinit(current_elem);
     112             : 
     113        1800 :     for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
     114             :     {
     115        1600 :       Real diff = std::abs(JxW_displaced[qp] - _JxW[qp]) / _JxW[qp];
     116        1600 :       if (diff > max_difference)
     117         581 :         max_difference = diff;
     118             :     }
     119             : 
     120         200 :     JxW_displaced.release();
     121         284 :   }
     122             : 
     123         284 :   _communicator.max(max_difference);
     124             : 
     125         284 :   if (max_difference > _max_jacobian_diff)
     126         189 :     return _max_jacobian_diff / max_difference;
     127             : 
     128             :   return 1.0;
     129         284 : }