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