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 "ADMortarLagrangeConstraint.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseVariable.h"
      14             : #include "MoosePreconditioner.h"
      15             : #include "VariableCondensationPreconditioner.h"
      16             : #include "NonlinearSystemBase.h"
      17             : #include "Assembly.h"
      18             : #include "SystemBase.h"
      19             : #include "ADUtils.h"
      20             : #include "FEProblemBase.h"
      21             : #include "AutomaticMortarGeneration.h"
      22             : 
      23             : #include <algorithm>
      24             : #include <vector>
      25             : 
      26             : InputParameters
      27        7465 : ADMortarLagrangeConstraint::validParams()
      28             : {
      29        7465 :   InputParameters params = ADMortarConstraint::validParams();
      30       14930 :   params.addParam<Real>(
      31             :       "derivative_threshold",
      32       14930 :       1.0e-17,
      33             :       "Threshold to discard automatic differentiation derivatives. This number is chosen on the "
      34             :       "order of the machine epsilon based on current experience.");
      35             : 
      36        7465 :   return params;
      37           0 : }
      38             : 
      39        4130 : ADMortarLagrangeConstraint::ADMortarLagrangeConstraint(const InputParameters & parameters)
      40             :   : ADMortarConstraint(parameters),
      41        8260 :     _ad_derivative_threshold(getParam<Real>("derivative_threshold")),
      42        4130 :     _apply_derivative_threshold(true)
      43             : {
      44        4130 : }
      45             : 
      46             : void
      47        4086 : ADMortarLagrangeConstraint::initialSetup()
      48             : {
      49        4086 :   SetupInterface::initialSetup();
      50             : 
      51             :   // Detect if preconditioner is VCP. If so, disable automatic derivative trimming.
      52        4086 :   auto const * mpc = feProblem().getNonlinearSystemBase(_sys.number()).getPreconditioner();
      53             : 
      54        4086 :   if (dynamic_cast<const VariableCondensationPreconditioner *>(mpc))
      55         118 :     _apply_derivative_threshold = false;
      56        4086 : }
      57             : 
      58             : void
      59    16775786 : ADMortarLagrangeConstraint::computeResidual(Moose::MortarType mortar_type)
      60             : {
      61    16775786 :   _primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
      62    16775786 :       *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
      63             : 
      64    16775786 :   _secondary_ip_lowerd_map = amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
      65             : 
      66             :   unsigned int test_space_size = 0;
      67    16775786 :   switch (mortar_type)
      68             :   {
      69     8384643 :     case Moose::MortarType::Secondary:
      70     8384643 :       prepareVectorTag(_assembly, _secondary_var.number());
      71     8384643 :       test_space_size = _test_secondary.size();
      72     8384643 :       break;
      73             : 
      74     8384643 :     case Moose::MortarType::Primary:
      75     8384643 :       prepareVectorTagNeighbor(_assembly, _primary_var.number());
      76     8384643 :       test_space_size = _test_primary.size();
      77     8384643 :       break;
      78             : 
      79        6500 :     case Moose::MortarType::Lower:
      80             :       mooseAssert(_var, "LM variable is null");
      81        6500 :       prepareVectorTagLower(_assembly, _var->number());
      82        6500 :       test_space_size = _test.size();
      83        6500 :       break;
      84             :   }
      85             : 
      86             :   std::vector<unsigned int> is_index_on_lower_dimension;
      87             :   // Find out if nodes are on the surface
      88   135704550 :   for (_i = 0; _i < test_space_size; _i++)
      89             :   {
      90   118928764 :     if (mortar_type == Moose::MortarType::Primary && !_primary_ip_lowerd_map.count(_i))
      91    29766012 :       continue;
      92             : 
      93    89162752 :     if (mortar_type == Moose::MortarType::Secondary && !_secondary_ip_lowerd_map.count(_i))
      94    29766012 :       continue;
      95             : 
      96    59396740 :     is_index_on_lower_dimension.push_back(_i);
      97             :   }
      98             : 
      99    76172526 :   for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
     100   281521464 :     for (const auto index : is_index_on_lower_dimension)
     101             :     {
     102   222124724 :       _i = index;
     103   444249448 :       _local_re(_i) += raw_value(_JxW_msm[_qp] * _coord[_qp] * computeQpResidual(mortar_type));
     104             :     }
     105             : 
     106    16775786 :   accumulateTaggedLocalResidual();
     107    16775786 : }
     108             : 
     109             : void
     110     5478274 : ADMortarLagrangeConstraint::computeJacobian(Moose::MortarType mortar_type)
     111             : {
     112             :   std::vector<ADReal> residuals;
     113             :   size_t test_space_size = 0;
     114             :   typedef Moose::ConstraintJacobianType JType;
     115             :   typedef Moose::MortarType MType;
     116             :   std::vector<JType> jacobian_types;
     117             :   std::vector<dof_id_type> dof_indices;
     118             :   Real scaling_factor = 1;
     119             : 
     120     5478274 :   _primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
     121     5478274 :       *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
     122             : 
     123     5478274 :   _secondary_ip_lowerd_map = amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
     124             : 
     125     5478274 :   switch (mortar_type)
     126             :   {
     127     2736487 :     case MType::Secondary:
     128     2736487 :       dof_indices = _secondary_var.dofIndices();
     129     2736487 :       jacobian_types = {JType::SecondarySecondary, JType::SecondaryPrimary, JType::SecondaryLower};
     130     2736487 :       scaling_factor = _secondary_var.scalingFactor();
     131     2736487 :       break;
     132             : 
     133     2736487 :     case MType::Primary:
     134     2736487 :       dof_indices = _primary_var.dofIndicesNeighbor();
     135     2736487 :       jacobian_types = {JType::PrimarySecondary, JType::PrimaryPrimary, JType::PrimaryLower};
     136     2736487 :       scaling_factor = _primary_var.scalingFactor();
     137     2736487 :       break;
     138             : 
     139        5300 :     case MType::Lower:
     140             :       mooseAssert(_var, "This should be non-null");
     141        5300 :       dof_indices = _var->dofIndicesLower();
     142        5300 :       jacobian_types = {JType::LowerSecondary, JType::LowerPrimary, JType::LowerLower};
     143        5300 :       scaling_factor = _var->scalingFactor();
     144        5300 :       break;
     145             :   }
     146             :   test_space_size = dof_indices.size();
     147             : 
     148     5478274 :   residuals.resize(test_space_size, 0);
     149             : 
     150             :   std::vector<unsigned int> is_index_on_lower_dimension;
     151             :   std::vector<dof_id_type> dof_indices_lower;
     152             : 
     153             :   // Find out if nodes are on the surface
     154             :   unsigned int number_indices_on_lowerd = 0;
     155             : 
     156    43291226 :   for (_i = 0; _i < test_space_size; _i++)
     157             :   {
     158    37812952 :     if (mortar_type == Moose::MortarType::Primary && !_primary_ip_lowerd_map.count(_i))
     159     9459216 :       continue;
     160             : 
     161    28353736 :     if (mortar_type == Moose::MortarType::Secondary && !_secondary_ip_lowerd_map.count(_i))
     162     9459216 :       continue;
     163             : 
     164    18894520 :     is_index_on_lower_dimension.push_back(_i);
     165    18894520 :     dof_indices_lower.push_back(dof_indices[_i]);
     166    18894520 :     number_indices_on_lowerd++;
     167             :   }
     168             : 
     169             :   std::vector<ADReal> residuals_lower;
     170     5478274 :   residuals_lower.resize(number_indices_on_lowerd, 0);
     171             : 
     172             :   // Only populate nodal residuals on the primary/secondary surfaces
     173             :   // We do this regardless of whether we are interpolating normals. Use of this class
     174             :   // implies we have Lagrange elements, so internal (high-dimensional) normals have no meaning
     175             :   // and should be zero. As such, we decide to omit them and avoid possible spurious population of
     176             :   // automatic differentiation-generated derivatives.
     177    24372794 :   for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
     178             :   {
     179             :     unsigned int index_lower = 0;
     180    88423944 :     for (const auto index : is_index_on_lower_dimension)
     181             :     {
     182    69529424 :       _i = index;
     183   139058848 :       residuals_lower[index_lower] += _JxW_msm[_qp] * _coord[_qp] * computeQpResidual(mortar_type);
     184             : 
     185             :       // Get rid of derivatives that we assume won't count (tolerance prescribed by user)
     186             :       // This can cause zero diagonal terms with the variable condensation preconditioner when the
     187             :       // no adaptivity option is used (dofs are not checked).
     188             :       // Uncomment when https://github.com/libMesh/MetaPhysicL/pull/18 makes it to MOOSE
     189             : 
     190    69529424 :       index_lower++;
     191             :     }
     192             :   }
     193             : 
     194     5478274 :   addJacobian(_assembly, residuals_lower, dof_indices_lower, scaling_factor);
     195     5478274 : }