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          38 : ElementJacobianDamper::validParams()
      22             : {
      23          38 :   InputParameters params = GeneralDamper::validParams();
      24          38 :   params.addClassDescription("Damper that limits the change in element Jacobians");
      25          76 :   params.addParam<std::vector<VariableName>>(
      26             :       "displacements", {}, "The nonlinear displacement variables");
      27          76 :   params.addParam<Real>(
      28             :       "max_increment",
      29          76 :       0.1,
      30             :       "The maximum permissible relative increment in the Jacobian per Newton iteration");
      31          38 :   params.set<bool>("use_displaced_mesh") = true;
      32          38 :   return params;
      33           0 : }
      34             : 
      35          19 : ElementJacobianDamper::ElementJacobianDamper(const InputParameters & parameters)
      36             :   : GeneralDamper(parameters),
      37          19 :     _tid(parameters.get<THREAD_ID>("_tid")),
      38          19 :     _assembly(_subproblem.assembly(_tid, _sys.number())),
      39          19 :     _qrule(_assembly.qRule()),
      40          19 :     _JxW(_assembly.JxW()),
      41          19 :     _fe_problem(*parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
      42          19 :     _displaced_problem(_fe_problem.getDisplacedProblem()),
      43          38 :     _max_jacobian_diff(parameters.get<Real>("max_increment"))
      44             : {
      45          19 :   if (_displaced_problem == NULL)
      46           0 :     mooseError("ElementJacobianDamper: Must use displaced problem");
      47             : 
      48          19 :   _mesh = &_displaced_problem->mesh();
      49             : 
      50          38 :   const std::vector<VariableName> & nl_vnames(getParam<std::vector<VariableName>>("displacements"));
      51          19 :   _ndisp = nl_vnames.size();
      52             : 
      53          76 :   for (unsigned int i = 0; i < _ndisp; ++i)
      54             :   {
      55          57 :     _disp_var.push_back(&_sys.getFieldVariable<Real>(_tid, nl_vnames[i]));
      56          57 :     _disp_incr.push_back(_disp_var.back()->increment());
      57             :   }
      58          19 : }
      59             : 
      60             : void
      61          19 : ElementJacobianDamper::initialSetup()
      62             : {
      63          19 : }
      64             : 
      65             : Real
      66         404 : ElementJacobianDamper::computeDamping(const NumericVector<Number> & /* solution */,
      67             :                                       const NumericVector<Number> & update)
      68             : {
      69             :   // Maximum difference in the Jacobian for this Newton iteration
      70         404 :   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             :   PARALLEL_TRY
      77             :   {
      78             :     try
      79             :     {
      80             :       // Loop over elements in the mesh
      81       20472 :       for (auto & current_elem : _mesh->getMesh().active_local_element_ptr_range())
      82             :       {
      83        9836 :         point_copies.clear();
      84        9836 :         point_copies.reserve(current_elem->n_nodes());
      85             : 
      86             :         // Displace nodes with current Newton increment
      87       88524 :         for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
      88             :         {
      89       78688 :           Node & displaced_node = current_elem->node_ref(i);
      90             : 
      91       78688 :           point_copies.push_back(displaced_node);
      92             : 
      93      314752 :           for (unsigned int j = 0; j < _ndisp; ++j)
      94             :           {
      95             :             dof_id_type disp_dof_num =
      96      236064 :                 displaced_node.dof_number(_sys.number(), _disp_var[j]->number(), 0);
      97      236064 :             displaced_node(j) += update(disp_dof_num);
      98             :           }
      99             :         }
     100             : 
     101             :         // Reinit element to compute Jacobian of displaced element
     102        9836 :         _assembly.reinit(current_elem);
     103        9832 :         JxW_displaced = _JxW;
     104             : 
     105             :         // Un-displace nodes
     106       88488 :         for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
     107             :         {
     108       78656 :           Node & displaced_node = current_elem->node_ref(i);
     109             : 
     110      314624 :           for (unsigned int j = 0; j < _ndisp; ++j)
     111      235968 :             displaced_node(j) = point_copies[i](j);
     112             :         }
     113             : 
     114             :         // Reinit element to compute Jacobian before displacement
     115        9832 :         _assembly.reinit(current_elem);
     116             : 
     117       88488 :         for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
     118             :         {
     119       78656 :           Real diff = std::abs(JxW_displaced[qp] - _JxW[qp]) / _JxW[qp];
     120       78656 :           if (diff > max_difference)
     121        1065 :             max_difference = diff;
     122             :         }
     123             : 
     124        9832 :         JxW_displaced.release();
     125         404 :       }
     126             :     }
     127           4 :     catch (MooseException & e)
     128             :     {
     129           0 :       _fe_problem.setException(e.what());
     130           0 :     }
     131           4 :     catch (std::exception & e)
     132             :     {
     133             :       // Allow the libmesh error/exception on negative jacobian
     134           4 :       const std::string & message = e.what();
     135           4 :       if (message.find("Jacobian") == std::string::npos)
     136           0 :         throw e;
     137             :       else
     138           4 :         _fe_problem.setException(message);
     139           4 :     }
     140             :   }
     141         404 :   PARALLEL_CATCH;
     142             : 
     143         398 :   _communicator.max(max_difference);
     144             : 
     145         398 :   if (max_difference > _max_jacobian_diff)
     146         297 :     return _max_jacobian_diff / max_difference;
     147             : 
     148             :   return 1.0;
     149         404 : }