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 "ReferenceElementJacobianDamper.h"
      11             : #include "FEProblem.h"
      12             : #include "MooseMesh.h"
      13             : #include "Assembly.h"
      14             : #include "SystemBase.h"
      15             : 
      16             : #include "libmesh/quadrature.h" // _qrule
      17             : 
      18             : registerMooseObject("SolidMechanicsApp", ReferenceElementJacobianDamper);
      19             : 
      20             : InputParameters
      21          28 : ReferenceElementJacobianDamper::validParams()
      22             : {
      23          28 :   InputParameters params = GeneralDamper::validParams();
      24          28 :   params += BlockRestrictable::validParams();
      25          28 :   params.addClassDescription("Damper that limits the change in element Jacobians");
      26          56 :   params.addRequiredCoupledVar("displacements", "The nonlinear displacement variables");
      27          56 :   params.addParam<Real>(
      28          56 :       "max_increment", 0.1, "The maximum permissible incremental Jacobian per Newton iteration");
      29          28 :   return params;
      30           0 : }
      31             : 
      32          14 : ReferenceElementJacobianDamper::ReferenceElementJacobianDamper(const InputParameters & parameters)
      33             :   : GeneralDamper(parameters),
      34             :     Coupleable(this, /*nodal=*/true),
      35             :     BlockRestrictable(this),
      36          14 :     _max_jacobian_diff(getParam<Real>("max_increment")),
      37          28 :     _tid(getParam<THREAD_ID>("_tid")),
      38          14 :     _mesh(_subproblem.mesh()),
      39          14 :     _assembly(_subproblem.assembly(_tid, _sys.number())),
      40          14 :     _qrule(_assembly.qRule()),
      41          14 :     _ndisp(coupledComponents("displacements")),
      42          42 :     _disp_num(coupledIndices("displacements"))
      43             : {
      44          14 :   _grad_phi.resize(_ndisp);
      45          56 :   for (auto i : make_range(_ndisp))
      46             :   {
      47          42 :     const auto & disp_var = _sys.getFieldVariable<Real>(_tid, _disp_num[i]);
      48          42 :     checkVariable(disp_var);
      49          42 :     _grad_phi[i] = &disp_var.gradPhi();
      50             :   }
      51          14 : }
      52             : 
      53             : Real
      54         240 : ReferenceElementJacobianDamper::computeDamping(const NumericVector<Number> & solution,
      55             :                                                const NumericVector<Number> & update)
      56             : {
      57             :   // Maximum difference in the Jacobian for this Newton iteration
      58         240 :   Real max_difference = 0.0;
      59             : 
      60             :   // Loop over elements in the mesh
      61        1008 :   for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
      62             :   {
      63         264 :     if (!hasBlocks(elem->subdomain_id()))
      64          88 :       continue;
      65             : 
      66             :     // Compute gradients of displacements before and after this update
      67         176 :     computeGradDisp(elem, solution, update);
      68             : 
      69             :     // Compute the element Jacobian difference
      70        1584 :     for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
      71             :     {
      72        1408 :       RankTwoTensor F;
      73        5632 :       for (auto i : make_range(_ndisp))
      74        4224 :         F.fillColumn(i, _grad_disp[i][qp]);
      75        1408 :       F.addIa(1);
      76             : 
      77        1408 :       RankTwoTensor F_update;
      78        5632 :       for (auto i : make_range(_ndisp))
      79        4224 :         F_update.fillColumn(i, _grad_disp_update[i][qp]);
      80        1408 :       F_update.addIa(1);
      81             : 
      82        1408 :       Real diff = std::abs(F_update.det() - F.det()) / F.det();
      83        1408 :       if (diff > max_difference)
      84         176 :         max_difference = diff;
      85             :     }
      86         240 :   }
      87             : 
      88         240 :   _communicator.max(max_difference);
      89             : 
      90         240 :   if (max_difference > _max_jacobian_diff)
      91         210 :     return _max_jacobian_diff / max_difference;
      92             : 
      93             :   return 1.0;
      94             : }
      95             : 
      96             : void
      97         176 : ReferenceElementJacobianDamper::computeGradDisp(const Elem * elem,
      98             :                                                 const NumericVector<Number> & solution,
      99             :                                                 const NumericVector<Number> & update)
     100             : {
     101             :   // Reinit variable shape functions
     102         176 :   _assembly.setCurrentSubdomainID(elem->subdomain_id());
     103         176 :   _assembly.reinit(elem);
     104             : 
     105         176 :   _grad_disp.resize(_ndisp);
     106         176 :   _grad_disp_update.resize(_ndisp);
     107         704 :   for (auto i : make_range(_ndisp))
     108             :   {
     109             :     std::vector<dof_id_type> dof_indices;
     110         528 :     _sys.dofMap().dof_indices(elem, dof_indices, _disp_num[i]);
     111         528 :     _grad_disp[i].resize(_qrule->n_points());
     112         528 :     _grad_disp_update[i].resize(_qrule->n_points());
     113        4752 :     for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
     114             :     {
     115        4224 :       _grad_disp[i][qp].zero();
     116             :       _grad_disp_update[i][qp].zero();
     117       38016 :       for (auto dof_idx : index_range(dof_indices))
     118             :       {
     119       33792 :         _grad_disp[i][qp] += (*_grad_phi[i])[dof_idx][qp] *
     120       33792 :                              (solution(dof_indices[dof_idx]) + update(dof_indices[dof_idx]));
     121       33792 :         _grad_disp_update[i][qp] += (*_grad_phi[i])[dof_idx][qp] * solution(dof_indices[dof_idx]);
     122             :       }
     123             :     }
     124         528 :   }
     125         176 : }