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         218 : LMKernel::validParams()
      21             : {
      22         218 :   auto params = ADKernelValue::validParams();
      23         436 :   params.addRequiredCoupledVar("lm_variable", "The lagrange multiplier variable");
      24         436 :   params.addParam<bool>(
      25             :       "lm_sign_positive",
      26         436 :       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         218 :   return params;
      31           0 : }
      32             : 
      33         121 : LMKernel::LMKernel(const InputParameters & parameters)
      34             :   : ADKernelValue(parameters),
      35         121 :     _lm_var(*getVar("lm_variable", 0)),
      36         121 :     _lm(adCoupledValue("lm_variable")),
      37         121 :     _lm_test(_lm_var.phi()),
      38         484 :     _lm_sign(getParam<bool>("lm_sign_positive") ? 1. : -1)
      39             : {
      40         121 : }
      41             : 
      42             : void
      43      249950 : LMKernel::computeResidual()
      44             : {
      45      249950 :   std::vector<Real> strong_residuals(_qrule->n_points());
      46             : 
      47      249950 :   precalculateResidual();
      48             : 
      49      999800 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      50      749850 :     strong_residuals[_qp] =
      51     1499700 :         MetaPhysicL::raw_value(this->precomputeQpResidual() * _ad_JxW[_qp] * _ad_coord[_qp]);
      52             : 
      53             :   // Primal residual
      54      249950 :   prepareVectorTag(_assembly, _var.number());
      55             : 
      56      999800 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      57     2999400 :     for (_i = 0; _i < _test.size(); _i++)
      58     2249550 :       _local_re(_i) += _test[_i][_qp] * strong_residuals[_qp];
      59             : 
      60      249950 :   accumulateTaggedLocalResidual();
      61             : 
      62             :   // LM residual
      63      249950 :   prepareVectorTag(_assembly, _lm_var.number());
      64             : 
      65      999800 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      66     2254650 :     for (_i = 0; _i < _lm_test.size(); _i++)
      67     1504800 :       _local_re(_i) += _lm_sign * _lm_test[_i][_qp] * strong_residuals[_qp];
      68             : 
      69      249950 :   accumulateTaggedLocalResidual();
      70      249950 : }
      71             : 
      72             : void
      73      150300 : LMKernel::computeResidualsForJacobian()
      74             : {
      75      150300 :   std::vector<ADReal> strong_residuals(_qrule->n_points());
      76      150300 :   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      150300 :   if (_residuals.size() != resid_size)
      81         121 :     _residuals.resize(resid_size, 0);
      82             : 
      83             :   // Resizing a std:vector doesn't zero everything
      84      903000 :   for (auto & r : _residuals)
      85      752700 :     r = 0;
      86             : 
      87      150300 :   _all_dof_indices = _var.dofIndices();
      88      150300 :   _all_dof_indices.insert(
      89      150300 :       _all_dof_indices.end(), _lm_var.dofIndices().begin(), _lm_var.dofIndices().end());
      90             : 
      91      150300 :   precalculateResidual();
      92             : 
      93      601200 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      94      901800 :     strong_residuals[_qp] = this->precomputeQpResidual() * _ad_JxW[_qp] * _ad_coord[_qp];
      95             : 
      96             :   // Primal
      97      601200 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      98     1803600 :     for (_i = 0; _i < _test.size(); _i++)
      99     2705400 :       _residuals[_i] += strong_residuals[_qp] * _test[_i][_qp];
     100             : 
     101             :   // LM
     102      601200 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     103     1356300 :     for (_i = 0; _i < _lm_test.size(); _i++)
     104     1810800 :       _residuals[_test.size() + _i] += _lm_sign * strong_residuals[_qp] * _lm_test[_i][_qp];
     105      150300 : }