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 "LMKernel.h"
      11             : #include "Assembly.h"
      12             : #include "MooseVariable.h"
      13             : #include "InputParameters.h"
      14             : #include "MooseArray.h"
      15             : #include "SystemBase.h"
      16             : 
      17             : #include "libmesh/quadrature.h"
      18             : 
      19             : InputParameters
      20         122 : LMKernel::validParams()
      21             : {
      22         122 :   auto params = ADKernelValue::validParams();
      23         244 :   params.addRequiredCoupledVar("lm_variable", "The lagrange multiplier variable");
      24         244 :   params.addParam<bool>(
      25             :       "lm_sign_positive",
      26         244 :       true,
      27             :       "Whether to use a positive sign when adding this object's residual to the Lagrange "
      28             :       "multiplier constraint equation. Positive or negative sign should be chosen such that the "
      29             :       "diagonals for the LM block of the matrix are positive");
      30         122 :   return params;
      31           0 : }
      32             : 
      33          65 : LMKernel::LMKernel(const InputParameters & parameters)
      34             :   : ADKernelValue(parameters),
      35          65 :     _lm_var(*getVar("lm_variable", 0)),
      36          65 :     _lm(adCoupledValue("lm_variable")),
      37          65 :     _lm_test(_lm_var.phi()),
      38         260 :     _lm_sign(getParam<bool>("lm_sign_positive") ? 1. : -1)
      39             : {
      40          65 : }
      41             : 
      42             : void
      43      168160 : LMKernel::computeResidual()
      44             : {
      45      168160 :   std::vector<Real> strong_residuals(_qrule->n_points());
      46             : 
      47      168160 :   precalculateResidual();
      48             : 
      49      672640 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      50      504480 :     strong_residuals[_qp] =
      51     1008960 :         MetaPhysicL::raw_value(this->precomputeQpResidual() * _ad_JxW[_qp] * _ad_coord[_qp]);
      52             : 
      53             :   // Primal residual
      54      168160 :   prepareVectorTag(_assembly, _var.number());
      55             : 
      56      672640 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      57     2017920 :     for (_i = 0; _i < _test.size(); _i++)
      58     1513440 :       _local_re(_i) += _test[_i][_qp] * strong_residuals[_qp];
      59             : 
      60      168160 :   accumulateTaggedLocalResidual();
      61             : 
      62             :   // LM residual
      63      168160 :   prepareVectorTag(_assembly, _lm_var.number());
      64             : 
      65      672640 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      66     1517520 :     for (_i = 0; _i < _lm_test.size(); _i++)
      67     1013040 :       _local_re(_i) += _lm_sign * _lm_test[_i][_qp] * strong_residuals[_qp];
      68             : 
      69      168160 :   accumulateTaggedLocalResidual();
      70      168160 : }
      71             : 
      72             : void
      73      101160 : LMKernel::computeResidualsForJacobian()
      74             : {
      75      101160 :   std::vector<ADReal> strong_residuals(_qrule->n_points());
      76      101160 :   const auto resid_size = _test.size() + _lm_test.size();
      77             :   mooseAssert(resid_size == (_var.dofIndices().size() + _lm_var.dofIndices().size()),
      78             :               "The test function sizes and the dof indices sizes should match.");
      79             : 
      80      101160 :   if (_residuals.size() != resid_size)
      81          65 :     _residuals.resize(resid_size, 0);
      82             : 
      83             :   // Resizing a std:vector doesn't zero everything
      84      607920 :   for (auto & r : _residuals)
      85      506760 :     r = 0;
      86             : 
      87      101160 :   _all_dof_indices = _var.dofIndices();
      88      101160 :   _all_dof_indices.insert(
      89      101160 :       _all_dof_indices.end(), _lm_var.dofIndices().begin(), _lm_var.dofIndices().end());
      90             : 
      91      101160 :   precalculateResidual();
      92             : 
      93      404640 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      94      606960 :     strong_residuals[_qp] = this->precomputeQpResidual() * _ad_JxW[_qp] * _ad_coord[_qp];
      95             : 
      96             :   // Primal
      97      404640 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      98     1213920 :     for (_i = 0; _i < _test.size(); _i++)
      99     1820880 :       _residuals[_i] += strong_residuals[_qp] * _test[_i][_qp];
     100             : 
     101             :   // LM
     102      404640 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     103      913320 :     for (_i = 0; _i < _lm_test.size(); _i++)
     104     1219680 :       _residuals[_test.size() + _i] += _lm_sign * strong_residuals[_qp] * _lm_test[_i][_qp];
     105      101160 : }