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