Line data    Source code

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