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        7878 : ADMortarLagrangeConstraint::validParams()
      28             : {
      29        7878 :   InputParameters params = ADMortarConstraint::validParams();
      30       15756 :   params.addParam<Real>(
      31             :       "derivative_threshold",
      32       15756 :       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        7878 :   return params;
      37           0 : }
      38             : 
      39        4345 : ADMortarLagrangeConstraint::ADMortarLagrangeConstraint(const InputParameters & parameters)
      40             :   : ADMortarConstraint(parameters),
      41        8690 :     _ad_derivative_threshold(getParam<Real>("derivative_threshold")),
      42        4345 :     _apply_derivative_threshold(true)
      43             : {
      44        4345 : }
      45             : 
      46             : void
      47        4299 : ADMortarLagrangeConstraint::initialSetup()
      48             : {
      49        4299 :   SetupInterface::initialSetup();
      50             : 
      51             :   // Detect if preconditioner is VCP. If so, disable automatic derivative trimming.
      52        4299 :   auto const * mpc = feProblem().getNonlinearSystemBase(_sys.number()).getPreconditioner();
      53             : 
      54        4299 :   if (dynamic_cast<const VariableCondensationPreconditioner *>(mpc))
      55         134 :     _apply_derivative_threshold = false;
      56        4299 : }
      57             : 
      58             : void
      59    21252138 : ADMortarLagrangeConstraint::computeResidual(Moose::MortarType mortar_type)
      60             : {
      61    21252138 :   _primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
      62    21252138 :       *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
      63             : 
      64    21252138 :   _secondary_ip_lowerd_map = amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
      65             : 
      66             :   unsigned int test_space_size = 0;
      67    21252138 :   switch (mortar_type)
      68             :   {
      69    10622049 :     case Moose::MortarType::Secondary:
      70    10622049 :       prepareVectorTag(_assembly, _secondary_var.number());
      71    10622049 :       test_space_size = _test_secondary.size();
      72    10622049 :       break;
      73             : 
      74    10622049 :     case Moose::MortarType::Primary:
      75    10622049 :       prepareVectorTagNeighbor(_assembly, _primary_var.number());
      76    10622049 :       test_space_size = _test_primary.size();
      77    10622049 :       break;
      78             : 
      79        8040 :     case Moose::MortarType::Lower:
      80             :       mooseAssert(_var, "LM variable is null");
      81        8040 :       prepareVectorTagLower(_assembly, _var->number());
      82        8040 :       test_space_size = _test.size();
      83        8040 :       break;
      84             :   }
      85             : 
      86    21252138 :   std::vector<unsigned int> is_index_on_lower_dimension;
      87             :   // Find out if nodes are on the surface
      88   174204186 :   for (_i = 0; _i < test_space_size; _i++)
      89             :   {
      90   152952048 :     if (mortar_type == Moose::MortarType::Primary && !_primary_ip_lowerd_map.count(_i))
      91    38277336 :       continue;
      92             : 
      93   114674712 :     if (mortar_type == Moose::MortarType::Secondary && !_secondary_ip_lowerd_map.count(_i))
      94    38277336 :       continue;
      95             : 
      96    76397376 :     is_index_on_lower_dimension.push_back(_i);
      97             :   }
      98             : 
      99    97649514 :   for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
     100   364706736 :     for (const auto index : is_index_on_lower_dimension)
     101             :     {
     102   288309360 :       _i = index;
     103   576618720 :       _local_re(_i) += raw_value(_JxW_msm[_qp] * _coord[_qp] * computeQpResidual(mortar_type));
     104             :     }
     105             : 
     106    21252138 :   accumulateTaggedLocalResidual();
     107    21252138 : }
     108             : 
     109             : void
     110     6748142 : 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     6748142 :   _primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
     121     6748142 :       *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
     122             : 
     123     6748142 :   _secondary_ip_lowerd_map = amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
     124             : 
     125     6748142 :   switch (mortar_type)
     126             :   {
     127     3370801 :     case MType::Secondary:
     128     3370801 :       dof_indices = _secondary_var.dofIndices();
     129     3370801 :       jacobian_types = {JType::SecondarySecondary, JType::SecondaryPrimary, JType::SecondaryLower};
     130     3370801 :       scaling_factor = _secondary_var.scalingFactor();
     131     3370801 :       break;
     132             : 
     133     3370801 :     case MType::Primary:
     134     3370801 :       dof_indices = _primary_var.dofIndicesNeighbor();
     135     3370801 :       jacobian_types = {JType::PrimarySecondary, JType::PrimaryPrimary, JType::PrimaryLower};
     136     3370801 :       scaling_factor = _primary_var.scalingFactor();
     137     3370801 :       break;
     138             : 
     139        6540 :     case MType::Lower:
     140             :       mooseAssert(_var, "This should be non-null");
     141        6540 :       dof_indices = _var->dofIndicesLower();
     142        6540 :       jacobian_types = {JType::LowerSecondary, JType::LowerPrimary, JType::LowerLower};
     143        6540 :       scaling_factor = _var->scalingFactor();
     144        6540 :       break;
     145             :   }
     146             :   test_space_size = dof_indices.size();
     147             : 
     148     6748142 :   residuals.resize(test_space_size, 0);
     149             : 
     150     6748142 :   std::vector<unsigned int> is_index_on_lower_dimension;
     151     6748142 :   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    54335474 :   for (_i = 0; _i < test_space_size; _i++)
     157             :   {
     158    47587332 :     if (mortar_type == Moose::MortarType::Primary && !_primary_ip_lowerd_map.count(_i))
     159    11903676 :       continue;
     160             : 
     161    35683656 :     if (mortar_type == Moose::MortarType::Secondary && !_secondary_ip_lowerd_map.count(_i))
     162    11903676 :       continue;
     163             : 
     164    23779980 :     is_index_on_lower_dimension.push_back(_i);
     165    23779980 :     dof_indices_lower.push_back(dof_indices[_i]);
     166    23779980 :     number_indices_on_lowerd++;
     167             :   }
     168             : 
     169     6748142 :   std::vector<ADReal> residuals_lower;
     170     6748142 :   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    30528122 :   for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
     178             :   {
     179             :     unsigned int index_lower = 0;
     180   112461240 :     for (const auto index : is_index_on_lower_dimension)
     181             :     {
     182    88681260 :       _i = index;
     183   177362520 :       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    88681260 :       index_lower++;
     191             :     }
     192             :   }
     193             : 
     194     6748142 :   addJacobian(_assembly, residuals_lower, dof_indices_lower, scaling_factor);
     195     6748142 : }