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