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 "SwitchingFunctionConstraintLagrange.h"
      11             : 
      12             : registerMooseObject("PhaseFieldApp", SwitchingFunctionConstraintLagrange);
      13             : 
      14             : InputParameters
      15         103 : SwitchingFunctionConstraintLagrange::validParams()
      16             : {
      17         103 :   InputParameters params = Kernel::validParams();
      18         103 :   params.addClassDescription("Lagrange multiplier kernel to constrain the sum of all switching "
      19             :                              "functions in a multiphase system. This kernel acts on the Lagrange "
      20             :                              "multiplier variable.");
      21         206 :   params.addParam<std::vector<MaterialPropertyName>>("h_names", "Switching function materials");
      22         206 :   params.addRequiredCoupledVar("etas", "eta order parameters");
      23         206 :   params.addParam<Real>("epsilon", 1e-9, "Shift factor to avoid a zero pivot");
      24         103 :   return params;
      25           0 : }
      26             : 
      27          55 : SwitchingFunctionConstraintLagrange::SwitchingFunctionConstraintLagrange(
      28          55 :     const InputParameters & parameters)
      29             :   : DerivativeMaterialInterface<JvarMapKernelInterface<Kernel>>(parameters),
      30         110 :     _h_names(getParam<std::vector<MaterialPropertyName>>("h_names")),
      31          55 :     _num_h(_h_names.size()),
      32          55 :     _h(_num_h),
      33          55 :     _dh(_num_h),
      34          55 :     _eta_map(getParameterJvarMap("etas")),
      35         165 :     _epsilon(getParam<Real>("epsilon"))
      36             : {
      37             :   // parameter check. We need exactly one eta per h
      38         110 :   if (_num_h != coupledComponents("etas"))
      39           0 :     paramError("etas", "Need to pass in as many etas as h_names");
      40             : 
      41             :   // fetch switching functions (for the residual) and h derivatives (for the Jacobian)
      42         197 :   for (std::size_t i = 0; i < _num_h; ++i)
      43             :   {
      44         142 :     _h[i] = &getMaterialPropertyByName<Real>(_h_names[i]);
      45             : 
      46         142 :     _dh[i].resize(_num_h);
      47         522 :     for (std::size_t j = 0; j < _num_h; ++j)
      48         760 :       _dh[i][j] = &getMaterialPropertyDerivative<Real>(_h_names[i], coupledName("etas", j));
      49             :   }
      50          55 : }
      51             : 
      52             : Real
      53     7789680 : SwitchingFunctionConstraintLagrange::computeQpResidual()
      54             : {
      55     7789680 :   Real g = -_epsilon * _u[_qp] - 1.0;
      56    29612240 :   for (std::size_t i = 0; i < _num_h; ++i)
      57    21822560 :     g += (*_h[i])[_qp];
      58             : 
      59     7789680 :   return _test[_i][_qp] * g;
      60             : }
      61             : 
      62             : Real
      63     4217520 : SwitchingFunctionConstraintLagrange::computeQpJacobian()
      64             : {
      65     4217520 :   return _test[_i][_qp] * -_epsilon * _phi[_j][_qp];
      66             : }
      67             : 
      68             : Real
      69    11251040 : SwitchingFunctionConstraintLagrange::computeQpOffDiagJacobian(unsigned int jvar)
      70             : {
      71    11251040 :   auto eta = mapJvarToCvar(jvar, _eta_map);
      72    11251040 :   if (eta >= 0)
      73             :   {
      74             :     Real g = 0.0;
      75    43481120 :     for (std::size_t i = 0; i < _num_h; ++i)
      76    32230080 :       g += (*_dh[i][eta])[_qp] * _phi[_j][_qp];
      77    11251040 :     return g * _test[_i][_qp];
      78             :   }
      79             : 
      80             :   return 0.0;
      81             : }